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Preface 


The  modeling  of  atmospheric  illumination  and  radiance  as  it  affects  Army 
systems  has  had  a  long  and  curious  history.  One  of  the  main  aspects  of  this  area 
has  been  the  lack  of  progress  in  solving  some  of  the  main  difficulties  in  assessing 
the  influence  of  clouds  and  hazes  in  complex  geometries.  The  example  of  the 
man  hunting  for  keys  beneath  the  street  lamp  comes  to  mind.  He  searched  where 
there  was  light,  but  his  keys  were  lost  someplace  else.  Often  in  the  atmospheric 
sciences  we  can  fall  victim  to  this  temptation  because  the  difficulty  of  some  real 
problems  is  beyond  current  theory  and  applications. 

In  this  case,  however,  we  believe  that  some  light  has  been  shed  on  the 
area  of  three-dimensional  modeling  of  radiances  in  cloud  fields.  And  while 
this  may  sound  esoteric,  the  issue  should  be  of  major  concern  to  the  Army 
community.  A  large  fraction  of  Army  systems  rely  on  direct  line  of  sight  (LOS) 
optics  for  engaging,  detecting,  and/or  tracking  enemy  activity.  And  since  the 
Army  frequently  operates  beneath  a  cloud  deck,  it  is  highly  likely  that  clouds 
will  impact  the  performance  of  Army  systems.  In  particular,  illumination 
effects  modify  the  ability  of  an  observer  to  detect  a  target  against  a  natural 
background.  Also,  in  the  past  the  illumination  and  optical  properties  along  lines 
of  sight  had  been  treated  as  uniform,  or,  in  the  best  models,  assumed  to  have 
simple  scaling  laws  for  use  in  slant-path  evaluation  techniques.  But,  in  three- 
dimensional  modeling,  the  radiance  held  changes  with  position  and  direction 
requiring  path  dependent  LOS  calculations.  This  allows  what-if  questions  to  be 
asked  concerning  into-sun/out-of-the-sun  scenarios  and  the  evaluation  of  cloud 
shadowing  effects.  Solar  illumination  also  affects  the  energy  loading  of  surfaces 
for  thermal  signature  calculations. 

The  Atmospheric  Illumination  Module  (AIM)  codes  are  designed  to  permit 
analysis  of  surface  illumination  and  volumetric  radiance  calculations  within  a 
scattering/absorbing/emitting  atmosphere  beginning  at  the  Earth’s  surface  and 
extending  upward  through  the  cloud-filled  portion  of  the  troposphere.  The 
principal  features  of  this  set  of  models  consist  of  a  three-dimensional  radiative 
transfer  model  capable  of  treating  thin  to  extremely  thick  media,  a  modeling 
methodology  for  representing  the  Legendre  expansion  of  the  aerosol  scattering 
properties  of  cloud  aerosols  and  hazes,  and  a  visualization  model  for  depicting 
the  appearances  of  clouds. 

It  is  hoped  the  visualization  model  discussed  herein  will  some  day  lead  to  the 
real  time  use  of  the  data  sets  produced  by  AIM  in  simulations  of  various  systems. 
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Executive  Summary 


Introduction 

The  Atmospheric  Illumination  Module  (AIM)  is  designed  to  facilitate  the 
assessment  of  radiance  in  the  atmosphere  (surface  through  4,  8,  or  12  km). 
Numerous  Army  systems  employ  passive  electro-optical  devices  designed  to 
detect  and  acquire  objects  in  their  image  fields,  both  those  directly  illuminated 
(visible)  and  those  detectable  through  indirect  effects  due  to  solar  loading.  But 
the  atmosphere  affects  the  appearance  of  viewed  scenes  both  through  spatial  and 
spectral  modifications  of  the  illumination  pattern  and  through  path  transmission 
and  radiance  modifications  of  the  scene  energy.  In  this  series  of  codes,  we 
currently  provide  a  means  of  visualizing  and  quantifying  the  effects  of  the 
atmosphere  at  visible  wavelengths. 


Background 

The  primary  impact  of  the  atmosphere  on  target  acquisition  in  the  visible  band 
is  contrast  transmission.  This  depends  on  both  the  density  of  the  medium 
(affecting  transmittance)  and  the  influence  of  illumination  on  atmospheric 
constituents  (affecting  path  radiance).  The  general  expression 


c 


Co 

l  +  5,(T-i-l) 


(1) 


is  often  used,  where  C  is  the  transmitted  contrast,  Cq  is  the  initial  contrast 
at  zero  range,  Sg  is  the  sky-to-ground  ratio,  and  T  is  the  transmission.  T  = 
exp(— cri?),  where  a  is  the  wavelength-dependent  extinction  coefficient  and  R 
is  the  range  to  target.  The  sky-to-ground  ratio  is  the  ratio  of  the  limiting 
path  radiance  (often  judged  based  on  the  brightness  of  the  horizon  sky)  to  the 
brightness  of  the  background  terrain  near  the  target. 

While  there  are  numerous  complications  that  apply  when  using  this  equation, 
as  discussed  in  the  main  text,  the  primary  point  is  that  limiting  path  radiance, 
terrain  brightness,  and  an  extinction  coefficient  are  necessary  for  any  analysis 
of  atmospheric  effects.  Providing  these  quantities  is  the  motivation  behind  the 
AIM  models. 
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Overview 


The  body  of  this  report  explains  the  methods  developed  for  producing  the 
volumetric  radiances  necessary  to  assess  sky  and  terrain  brightnesses.  We 
consider  the  means  of  generating  suitable  model  inputs  by  converting  liquid 
water  content  information  provided  by  the  Air  Force’s  Cloud  Scene  Simulation 
Model  (CSSM)  into  extinction  coefficient  and  scattering  property  data  required 
by  the  radiative  transfer  (RT)  model. 

Assessing  brightness  entails  the  RT  modeling  of  a  simulated  atmospheric  volume 
in  contact  with  the  Earth’s  surface  and  extending  several  kilometers  vertically 
to  incorporate  typical  cloud  layers  for  surface  shadow  calculation  purposes.  The 
scenario  volume  is  extended  horizontally  to  treat  a  volume  large  enough  to  be 
significant  for  line  of  sight  (LOS)  engagements  (around  8x8  km^). 

The  RT  model  developed  here  features  techniques  for  treating  dense  clouds 
within  a  Discrete  Ordinates  Method  (DOM)  paradigm.  This  is  accomplished  by 
altering  the  traditional  DOM  technique.  The  traditional  technique  emphasizes 
cell-centered  calculations  and  extrapolations  to  cell-wall  values.  The  modified 
technique  determines  cell-wall-averaged  radiances  and  extrapolates  results  to 
cell- volume-averaged  radiances  for  purposes  of  calculating  volumetric  diffuse 
scattering  effects.  After  these  volume-based  methods  were  developed,  analysis 
showed  that  energy  was  being  lost  due  to  the  single-scattering  assumption.  To 
compensate,  a  surface-based  scattering  event  has  been  introduced,  based  on  a 
volumetric  energy  accounting  calculation.  Validations  of  the  RT  code  have  been 
accomplished  through  comparison  with  Monte  Carlo  scenario  results.  These 
comparisons  involved  developing  a  series  of  three-dimensional  scattering  cases, 
which  in  one  limit  matched  existing  literature  results  for  uniform  density  cubical 
volumes. 

Following  development  of  the  RT  technique,  we  consider  the  usages  of  these 
results  in  LOS  effects  computations,  visualization  codes,  and  compact  formats 
useful  for  network  distribution  of  data.  In  the  appendix  we  show  the  means 
of  representing  aerosol  scattering  properties  within  the  Legendre  expansion 
framework  applicable  to  the  DOM. 


Conclusions 

Methods  used  in  developing  radiance  mappings  and  visualization  techniques  for 
the  near-surface  atmosphere  for  the  visible  band  are  described.  The  RT  model 
developed  is  shown  to  predict  model  fluxes  with  0.999  correlation  to  the  Monte 
Carlo  cases  run.  The  visualization  code  developed  is  shown  to  produce  realistic 
images  that  are  also  quantitatively  correct.  We  believe  these  models  represent 
the  state  of  the  art  in  physics-based  radiance  modeling  of  the  boundary  layer 
atmosphere. 
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1.  Introduction 


The  performance  of  any  electro-optical  imaging  system  depends  on  the  medium 
through  which  signals  are  received.  As  optical  systems  are  continually  being 
upgraded  for  special  purposes,  it  is  reasonable  to  augment  system  design  by 
studying  the  impact  of  environment  on  system  performance.  One  means  of 
studying  this  problem  is  to  simulate  the  system  along  with  its  environment 
in  such  a  way  that  system-controlling  parameters  can  be  modified  and  tested 
against  a  variety  of  adverse  situations.  The  optimal  performance  can  then  be 
determined  as  that  combination  of  parameters  providing  the  best  functionality 
under  the  conditions  tested. 

1.1  Target  Acquisition  Considerations 

From  a  systems  viewpoint,  the  primary  effect  of  the  atmosphere  is  the  loss  of 
contrast  between  a  target  and  its  background.  An  Army  adage  says,  “if  it  can 
be  seen  it  can  be  killed.”  But  what  does  it  mean  for  something  to  be  seen  using 
a  visible  waveband  system? 

Middleton  (1952)  presented  the  classic  study  indicating  the  human  vision  system 
has  a  limiting  contrast  threshold  of  approximately  2  percent  under  daylight 
conditions.  Contrast  is  usually  defined  as 


C{R) 


Lt{R)  -  Lb{R) 
Lb{R) 


(2) 


where  C  is  the  contrast  and  Lt  and  Lb  are  the  propagated  brightnesses 
(radiances)  of  the  target  (object)  and  its  background,  all  at  range  R.  When 
C{R)  falls  below  0.02,  the  object  becomes  indistinguishable  from  its  background 
to  a  human  observer. 

Standard  analysis  of  Eq.  (2)  involves  defining  Cq  as  the  contrast  at  zero  range 
and  a  as  the  wavelength-dependent  extinction  coefficient  (assumed  constant  over 
the  optical  path  to  the  observer).  The  initial  radiance  of  the  target  is  reduced  by 
scattering  and  absorption  along  the  optical  path  by  a  factor  T(R)  =  exp(—aR) 
called  the  transmittance.  In  addition,  a  quantity  (Lp(R))  known  as  the  path 
radiance  (scattering  and  emission  of  energy  into  the  LOS)  is  added  to  both  the 
attenuated  target  and  background  radiances: 


Lt{R)  =  T{R)Lt{Q)  +  Lp{R)-  Lb{R)  =  T{R)Lb{Q)  +  Lp{R).  (3) 
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In  the  traditional  analysis,  the  path  radiance  is  related  to  the  quantity  called 
the  limiting  path  radiance  or  sky  radiance  (-h^),  which  is  the  radiance  produced 
if  one  could  look  through  an  infinitely  long  atmosphere  in  a  given  direction.  The 
relationship  between  the  path  radiance  and  the  limiting  path  radiance  is: 


Lp{R)  =  Ll  [1-T{R)]. 


(4) 


Introducing  Eqs.  (2)  and  (3)  into  (1),  we  obtain  the  form, 


^  T{R)ILt{0)  -  Ib(0)] 

^  '  r(B)LB(0)  +  ii[i-r(ij)i' 


(5) 


Dividing  numerator  and  denominator  by  T(R)Lb{0)^  and  using  the  terms  sky- 
to-ground  ratio  (Sg  =  Ll / Lb{0))  and  inherent  contrast  (Co  =  [Lt{0)  — 
Lb{0)]/ Lb{0))^  one  obtains. 


C{R) 


C„ 

1  +  S,{T{R)-^ -ly 


(6) 


which  is  the  contrast  transmission  equation.  This  equation  has  important 
consequences  in  assessing  atmospheric  effects  on  target  acquisition  and  thus 
has  been  incorporated  into  numerous  combat  simulations.  The  influence  of 
transmittance  is  obvious.  The  influence  of  sky-to-ground  ratio  is  less  direct, 
but  not  negligible.  Lee  and  LaMotte  (1991)  have  shown  how  variations  in  Sg 
influence  the  loss  exchange  ratios  in  combat  simulations. 

But  there  are  several  difficulties  in  the  use  of  this  simple  expression  in  describing 
the  effects  of  the  atmosphere.  First,  Cq  is  by  no  means  ‘inherent’  in  the  sense 
that  it  is  an  invariant  property  of  a  given  target  against  a  given  background. 
The  appearance  of  an  object  against  its  background  will  depend  on  the  available 
direct  (solar/lunar)  and  diffuse  illumination’s  interactions  with  an  object  or 
background’s  bi-directional  reflectance  distribution  function  (BRDF)  (Davis 
(1987);  Weiss  and  Scoggins  (1987)),  or  even  its  facetized  representation  in  more 
complex  modeling  schemes.  These  characteristics  may  also  be  spectrally  (color) 
dependent  (e.  g.,  Gerhart  et  al.  (1995)). 

Second,  transmittance  is  not  a  simple  function  of  range.  Slant  path  geometries 
through  the  atmosphere  can  involve  vertically  variable  extinction  coefficients 
through  haze  (Fiegel  1994)  or  fully  three-dimensional  (3D)  variability  when 
considering  clouds. 

Third,  the  sky  brightness  or  limiting  path  radiance  is  a  function  of  not  only 
position  and  wavelength  but  also  of  look  direction  (Davis  1986)  due  to  non¬ 
isotropic  scattering.  Further,  dense  volume  segments  will  create  dark  and  bright 
shadow  regions  in  the  volume,  invalidating  the  simple  form  for  the  path  radiance 
calculation.  To  properly  evaluate  a  radiative  environment,  then,  one  needs  a 
fully  3D  spectral/directional  model. 
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1.2  The  Modeling  Environment 

These  same  observations  are  made  in  various  studies  (Welch  and  Wielicki  (1984); 
Kobayashi  (1991);  Haferman  et  al.  (1993);  Li  eA  al.  (1994);  Byrne  eA  al. 
(1996))  which  consider  the  ambient  radiances  present  within  the  tropospheric 
atmosphere  and  atmospheric  layer  reflectivities  of  partly  cloudy  atmospheres. 
These  studies  concluded  that  one  cannot  adequately  model  partly  cloudy  skies 
using  plane-parallel  atmospheric  models,  where  typically  partly  cloudy  results 
are  interpolated  from  results  of  overcast  and  clear-sky  cases.  These  studies 
thus  infer  the  implicit  weaknesses  of  plane-parallel  models  to  determine  target 
acquisition  capabilities  within  the  Earth’s  cloud-filled  atmosphere.  Yet  plane- 
parallel  models  still  have  numerous  applications.  For  example,  plane-parallel 
radiative  transfer  treatments  heavily  influence  considerations  of  vertical  remote 
sensing  through  cloud  free  atmospheres  (Nakajima  and  Tanaka  1988)  and  high 
accuracy  one-dimensional  (ID)  radiative  transfer  (RT)  codes  (Stamnes  and 
Swanson  (1981);  Stamnes  eA  al.  (1988))  have  been  developed  expressly  for 
analysis  of  these  problems. 

While  ID-RT  codes  are  clearly  valuable  for  certain  applications,  robust  analysis 
of  systems  performance  within  a  realistic  earth  atmosphere  requires  a  3D-RT 
code.  Of  the  numerous  3D-RT  techniques  described  in  the  literature,  several  are 
generally  recognized,  including  Fourier  methods,  analytical  models,  diffusion 
approximation  models,  Monte  Carlo  models,  and  Discrete  Ordinates  Methods 
(DOM)  models.  Fourier  method  studies  are  typified  by  Priesendorfer  and 
Stephens  (1984),  Stephens  (1986,  1988),  Kobayashi  (1991),  and  Li  et  al.  (1994b). 
Analytical  methods  have  been  developed  by  Davies  (1978),  used  by  van  de 
Hulst  (1980a,  1980b),  and  developed  for  special  cases  more  recently  by  Li  et  al. 
(1994a).  Diffusion  methods  were  pursued  by  Cube  et  al.  (1980)  and  Zardecki 
et  al.  (1986).  Monte  Carlo  methods  were  used  by  McKee  and  Cox  (1974) 
and  Welch  and  Wielicki  (1984).  Yet  another  approach  involves  the  small-angle- 
approximation  (SAA).  This  method  is  often  used  to  study  systems  effects  due  to 
forward  scattering  of  aerosols.  But,  we  are  considering  only  general  illumination 
and  visualization  of  scenes  and  the  resulting  dynamic  range  effects.  In  effect, 
we  are  ignoring  the  forward  scatter  problem,  which  is  a  separate  topic,  and  thus 
will  ignore  the  SAA  approach  here. 

For  this  work  we  have  selected  the  discrete-ordinates  technique  as  the  most 
appropriate  vehicle.  This  focus  is  prompted  by  several  factors.  First,  diffusion 
methods  are  inherently  limited  to  relatively  isotropically  scattering  media. 
Moreover,  this  media  must  be  optically  thick  and  have  extinction  effects  that 
are  dominated  by  scattering.  The  atmosphere  viewed  by  sensors  does  not  fit 
this  description  in  general;  while  diffusion  processes  may  dominate  in  extremely 
thick  regions  (clouds)  of  the  atmosphere,  they  have  no  applicability  to  non¬ 
cloud  regions.  But,  any  general  model  must  be  able  to  treat  both  thin  and  thick 
regions. 
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Monte  Carlo  techniques  are  currently  impractical,  or  at  best  unwieldy,  given 
current  limited  computer  resources.  In  order  to  characterize  paths  through 
the  volume,  one  must  produce  a  radiance  picture  of  the  whole  volume,  and 
Monte  Carlo  methods  do  not  work  well  at  characterizing  all  locations  in  all 
directions.  Similar  considerations  restrict  the  usage  of  Fourier  methods,  where 
large  matrices  must  be  inverted  to  produce  a  solution.  Typical  studies  such 
as  those  by  Kobayashi  (1991)  and  Li  et  al.  (1994b)  restrict  the  positional 
or  angular  variation  of  the  analysis  space  in  order  to  limit  the  size  of  the 
matrix  to  be  inverted.  For  the  size  of  our  sample  volume,  at  a  resolution 
that  characterizes  typical  natural  cloud  dimensions,  the  size  of  matrix  necessary 
becomes  impossibly  large.  For  example,  Fiveland  (1987)  indicates  at  least  a 
24-stream  solution  is  required  to  provide  reasonable  directional  information  on 
radiant  intensity  for  visualization  purposes.  To  achieve  this  resolution  for  our 
typical  32^-cell  spatial  model  using  a  Fourier  method,  one  would  need  to  invert 
a  matrix  of  (32^  X  24)^  =  786,432^  =  6.2  x  10^^  elements.  This  approach 
appears  highly  impractical.  Other  analytical  methods  appear  to  be  as  restricted 
in  application  as  the  Fourier  methods.  These  considerations  led  to  the  selection 
of  the  DOM  in  3D  using  a  finite-element  approach. 

DOM  was  originally  developed  by  Chandrasekhar  (1960),  from  which  two 
parallel  courses  were  pursued.  The  first,  primarily  useful  in  astronomical 
and  satellite  applications,  focused  on  ID  modeling  (Lion  (1973),  Stamnes  and 
Swanson  (1981),  Stamnes  eA  al.  (1988)),  culminating  in  the  development 
of  the  DISORT  (DIScrete  Ordinates  Radiative  Transfer)  routine,  a  standard 
in  the  field.  The  second  effort  focused  on  3D  modeling,  with  applications 
primarily  in  furnace  technology  and  nuclear  scattering.  The  basic  theory  for  this 
application  was  the  so-called  Finite  Element  (FE)  DOM  (Carlson  and  Lathrop 
1968).  Improvements  to  this  theory  included  better  simulation  of  scattering 
via  the  scaling  group  transformation  (McKellar  and  Box  1981),  as  commonly 
implemented  through  Wiscombe’s  (1977)  (5-M  method.  Numerous  applications 
of  the  DOM  include  study  of  radiative  interaction  with  the  Earth’s  surface 
(Zardecki  et  al.  (1983),  Gerstl  and  Zardecki  (1985a,  1985b)),  furnace  technology 
(Eiveland  (1985,  1987)),  and  foliated  canopies  (Myneni  et  al.  (1990,  1991)).  But, 
problems  still  remain  with  handling  dense  media  under  DOM.  Some  attempts 
have  been  made  to  correct  these  problems  (Wakil  and  Sakadura  (1992),  Chai  et 
al.  (1994)),  but  these  have  mainly  been  patches  to  the  basic  theory. 

The  development  described  here  has  its  roots  in  neutron  scattering  work  of  the 
70’s.  We  use  a  level  symmetric  even  quadrature  set  (Lathrop  and  Carlson  (1965); 
Eiveland  (1991))  proposed  by  Lewis  and  Miller  (1984).  More  recently,  Zardecki 
and  Davis  (1991)  and  Wetmore  and  Zardecki  (1993)  developed  the  immediate 
predecessor  of  the  model  described  herein:  the  Boundary  Layer  Illumination  and 
Radiative  Balance  model  (BLIRB). 

Unfortunately  BLIRB  attempted  to  model  3D  effects  using  traditional  EE  DOM, 
leading  to  the  use  of  the  so-called  negative  flux  fixup  patch.  The  problem 
of  negative  fluxes  (Eiveland  (1985);  Gerstl  and  Zardecki  (1985a);  Wakil  and 
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Sakadura  (1992)),  occurs  when  the  optical  thickness  across  model  cells  exceeds  a 
threshold  value.  Then,  the  computational  method  for  evaluating  stream  radiance 
exiting  the  cell  can  result  in  negative  values.  Under  the  traditional  approach, 
the  only  means  of  avoiding  this  problem  is  to  limit  the  maximum  optical  depth 
of  a  cell  by  limiting  the  cell  dimensions  or  through  use  of  an  auxiliary  method 
of  compensation  known  as  the  negative  flux  fixup.  This  coding  patch  artificially 
allows  denser  media  within  a  simulation,  but  it  does  not  alleviate  the  consequent 
lack  of  accuracy  of  the  model  results.  That  is,  how  do  you  establish  what  the 
right  value  should  be  if  your  estimate  is  a  negative  number? 

Although  alternative  corrections  have  been  suggested  (Wakil  and  Sacadura 
1992),  a  general  FE  DOM  for  arbitrary  media  density  has  heretofore  been 
lacking.  Moreover,  DOM  techniques  have  often  employed  flux  renormalizations 
to  conserve  energy  (Myneni  et  al.  1990).  But  since  the  negative  flux  fixup 
produces  cloud  regions  that  are  artificially  darkened,  renormalization  attempts 
then  merely  tend  to  overbrighten  regions  between  clouds  without  properly 
restoring  cloud  brightness. 

In  section  2,  a  methodology  for  alleviating  this  modeling  deficiency  is  presented 
where  a  surface-based  DOM  computation  technique  is  proposed.  In  section  3,  a 
Monte  Carlo  model  is  described  and  compared  with  standard  literature  results 
(McKee  and  Cox  1974).  In  section  4,  this  Monte  Carlo  model  is  used  to  evaluate 
the  results  of  the  proposed  method,  using  the  McKee  and  Cox  scenarios  and  3D 
extensions  to  these  scenarios.  In  section  5,  the  results  of  the  radiative  transfer 
code  are  utilized  in  a  cloud  rendering  technique,  and  the  formats  of  standard 
outputs  are  described  for  use  in  other  applications. 

1.3  Support 

This  software  was  developed  as  a  consequence  of  funding  under  several  joint 
programs  mentioned  in  the  acknowledgements  section.  Currently,  products  from 
this  software  are  provided  through  the  Master  Environmental  Library  project, 
a  joint  DoD  program  sponsored  by  the  Defense  Modeling  and  Simulation  Office. 
Support  for  this  software  can  be  obtained  by  contacting: 

1.3.1  Mailing  Address 

U.S.  Army  Research  Laboratory 
ATTN:  AMSRL-IS-EW  (D.  Tofsted) 

White  Sands  Missile  Range,  NM  88002-5501 
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1.3.2  Phone  Numbers  and  E-Mail  Addresses 


Radiative  transfer  models  and  derived  output  products: 

COM:(505)  678-3039 
FAX:(505)  678-3385 
DSN:  258-3039 

emaihdtof  stedOarl  .mil 

Visualization  products  and  LOS  integrations  using  RT  code  outputs: 
COM:(505)  678-1570 
FAX:(505)  678-3385 
DSN:  258-1570 

emaiLsobrienOpsl .  nmsu .  edu 
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2.  DOMRT 


2.1  The  Radiative  Transfer  Equation 

The  Discrete-Ordinate  Method  approach  to  Radiative  Transfer  (DOMRT)  is  a 
standard  technique  used  in  numerically  solving  the  equation  of  radiative  transfer. 
We  begin  by  describing  the  traditional  approach. 

In  the  initial  step,  the  scattering  volume  is  divided  into  a  collection  of  gridded 
cells.  Each  cell  is  assumed  to  exhibit  uniform  scattering  properties  throughout. 
The  method  then  determines  the  radiance  properties  of  diffuse  radiation  flowing 
within  each  cell  in  a  series  of  directions.  In  deriving  the  standard  equation  used, 
one  begins  with  the  equation  of  radiative  transfer. 


O  •  V/(F,  O)  -|-  (7  /(F,  O)  =  as 


I{  f,  O')  P(0,  O')  dO'  +  (a  -  as)h{T{  f),  A),  (7) 

47r 


where,  /(F,  O)  is  the  streaming  radiance  at  position  F  and  in  direction  O.  O  is 
a  unit  vector  composed  of  direction  cosines  /i,  ?/,  and  along  the  x,  y,  and  2; 
axes,  respectively.  P(0,  O')  is  referred  to  as  the  phase  function,  representing  the 
probability  of  scattering  of  a  photon  between  incident  and  departing  directions 
(O'  and  O,  respectively),  given  that  a  scattering  occurs.  The  normalization 
condition  for  P(0,0')  is  that  integration  over  all  possible  scattering  directions 
equals  unity: 

/  P(0,  0')d0' =  1.  (8) 

J  47T 

Also,  due  to  the  principle  of  reciprocity  (Van  de  Hulst  1980a),  one  must  have 

P(0,  O')  =  P(0',  O).  (9) 

The  first  term  on  the  left  in  Eq.  (7)  represents  the  differential  change  in  the 
streaming  radiance,  I,  along  direction  O.  The  second  term  on  the  left  represents 
combined  extinction  losses  in  the  stream  radiance  due  to  absorption  and 
scattering.  The  right-hand  side  (RHS)  consists  of  contributions  to  the  stream 
radiance  due  to  scattering  and  graybody  emissions  within  the  medium  (first 
and  second  terms,  respectively).  In  the  scattering  integral,  dO'  is  a  differential 
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solid  angle,  a  is  the  extinction  coefficient,  ag  is  the  scattering  coefficient,  and 
b{T{r),  A)  is  the  blackbody  emission  function. 

Let  us  now  simplify  the  writing  of  Eq.  (7)  through  the  use  of  the  terms  J  and 
i?,  where  J  is  the  scattering  source  term. 


J(0,  r)  =  w{f)  /  /(r.  O')  P(0,  O')  dO', 

J  47T 

and  B  is  the  emissive  source  term, 

B{r,  A)  =  (1  -  w{r))  b{K{f),  A), 


(10) 


(11) 


where  A  is  the  radiation  wavelength  (expressed  here  in  /xm),  and  K  is  the 
temperature  (in  Kelvin)  of  the  medium,  as  a  function  of  position.  The  coefficient 
ru  is  the  single  scattering  albedo,  representing  the  probability  that  a  collision 
of  a  photon  with  a  particle  will  result  in  a  scattering  event.  The  assumption 
of  local  thermodynamic  equilibrium  dictates  that  the  probability  of  an  emission 
event  is  proportional  to  {1  —  zu).  Hence,  a  —  as  =  {l  —  zu)a  =  a  a  is  the  absorption 
coefficient.  Together,  the  sum  of  J  and  B  equals  the  limiting  path  radiance,  L. 


The  blackbody  emission  function  may  be  written  as, 

1.19106  xlO«A-5 

6(T(r),  A)  = 


exp[14,388/(AT)]  -  1’ 


(12) 


producing  a  result  with  units  of  W/m^-sr-/im.  (All  results  of  this  derivation 
represent  spectral  radiances,  but  the  frequency  dependence  hereafter  will  be 
suppressed.) 

Using  these  definitions,  the  radiative  transfer  equation  reduces  to. 


O  •  V/(F,  O)  +  (7  /(F,  O)  =  a  J(  f^  O)  +  B(r) 


(13) 


If  the  medium  is  now  discretized  such  that  T,  cr,  ru,  and  the  scattering  source 
integral  are  assumed  constant  over  the  span  of  a  cell,  then  the  RHS  of  Eq.  (13) 
is  constant,  and  one  can  introduce  the  variable  T,  defined  as, 

T(0,  r)  =  /(O,  r)  -  T(0,  r)  =  /(O,  r)  -  J(0,  r)  -  H(r,  A).  (14) 

Introducing  this  variable  in  Eq.  (13)  we  find. 


O- VT(r,0)  +  (7T(r,0)  =  0.  (15) 

This  equation  can  be  solved  if  we  make  two  more  assumptions.  Eirst,  we  assume 
that  the  input  stream  values  to  a  cell  are  constant  over  each  input  face.  Second, 
we  evaluate  the  scattering  source  term  J  by  taking  an  average  of  the  unscattered 
stream  energy  over  the  volume  of  each  cell.  We  then  evaluate  the  average 
radiance  of  a  given  stream  exiting  a  given  cell  face  by  averaging  the  output 
results  over  that  face.  This  average  then  forms  the  basis  for  use  as  an  incident 
value  in  the  next  adjacent  cell. 
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2.2  Level  Symmetry  Quadrature  Sets 

Before  continuing  with  this  train  of  thought,  though,  we  must  consider  the 
heart  of  the  discrete  ordinates  approach.  Under  this  approach,  all  integrations 
over  solid  angle  are  approximated  by  summations  over  a  specially  selected  set 
of  sample  directions,  (Sets  of  sample  directions  are  available  at  increasing 

numbers  of  sample  stream  points.)  As  higher  order  ordinate  expansions  are  used, 
the  method  becomes  progressively  better.  In  our  study,  level  symmetric  sets  of 
ordinates  using  Gaussian  quadrature  have  been  chosen  (Lewis  and  Miller  1984). 
While  these  sets  have  certain  limitations  regarding  flux  integration  (Fiveland 
1991),  we  minimize  the  error  associated  with  these  integrations  in  the  approach 
developed. 

A  level  symmetry  set  means  that  the  diffuse  stream  directions  are 

characterized  by  combinations  of  a  single  set  of  component  values  Xi  which 
are  symmetric  about  the  zero  planes  in  x,  y,  and  2;. 

Each  order  of  approximation  is  denoted  by  the  symbol  Sjv,  where  S  indicates  a 
level  symmetry  and  N  represents  the  order  of  approximation.  The  values  for  N 
are  even  integers  ranging  upward  from  2.  For  a  given  order  of  approximation  iV, 
there  will  be  iV  X  {N  +  2)  =  M  streams.  Thus  S2  is  an  8-stream  approximation, 
S4  is  a  24-stream  model,  and  Se  is  a  48-stream  model.  For  the  time  being  we 
will  focus  on  the  S2  and  S4  cases,  though  extensions  to  higher  orders  are  readily 
available. 

To  understand  how  the  number  of  streams  are  related  to  the  order  of 
approximation,  consider  that  each  streaming  direction  can  be  described 
by  a  set  of  three  integers  z,  j.  A;,  where  y  =  yp  y  =  Xj,  C  =  X/t,  such  that 

X?  +  Xj  +  xl  =  1-  (16) 

There  are  iV/2  —  1  unique  triplets  for  level  symmetric  set  Sjv-  Lathrop  and 
Carlson  (1965)  have  shown  that  only  one  degree  of  freedom  remains  in  choosing 
the  direction  cosines,  using  the  relation: 

X  =  \;+2^^(  1-3x5).  (17) 

Of  course,  when  N  =  2  there  are  no  degrees  of  freedom,  since,  necessarily, 
Xi  =  UTS;  X2  =  ~Xi-  III  llic  derivation  of  Eq.  (17)  the  last  term  is  assumed 
positive,  which  requires  0  <  Xi  A  1/3  (Fiveland  1991). 

For  the  S2  and  S4  approximations  considered  in  this  paper,  the  Xi’s  are  —^1/3 
and  x/T/S  in  the  8-stream  approximation,  and  equal  to  -0.8688903,  -0.3500212, 
0.3500212,  and  0.8688903  in  the  24-stream  model.  Thus,  for  S2,  there  is  only 
one  combination  of  x  values  possible  for  an  O  vector  in  the  first  octant,  or 
eight  possible  directions  when  all  octants  are  considered.  For  S4,  there  are 
three  possibilities  for  vectors  in  the  first  octant:  let  xi  =  0.3500212  and  X2  = 
0.8688903.  Then,  Oi  =  (x2,  Xi,  Xi);  ^2  =  (xi,  X2,  Xi);  ^3  =  (xi,  Xi,  X2). 
Similar  permutations  occur  in  each  of  the  other  seven  octants,  resulting  in  24 
streams  total. 
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2.3  Angular  Integration 

The  value  of  the  discrete  ordinates  approach  is  that  scattered  energy  at  any 
angle  can  be  interpolated  to  the  degree  of  resolution  of  the  approximation  used. 
Thus,  there  is  always  a  means  of  expressing  radiance  in  all  directions,  as  opposed 
to  methods  such  as  Monte  Carlo,  where  the  means  of  interpolating  intermediate 
directions  can  be  ambiguous. 

The  general  equation  for  performing  angular  integrations  under  the  discrete 
ordinates  approach  is  via  the  approximation: 
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where  are  the  discrete  ordinate  directions,  which  depend  on  the  particular 
order  of  approximation  and  ordinate  direction  selection  scheme  chosen,  and  the 
Wm  are  weighting  coefficients.  In  essence,  the  measure  the  fraction  of  an 
octant  of  solid  angle  that  is  being  occupied  (influenced)  by  that  ordinate;  in 
particular,  =  1  in  the  8-stream  model  and  =  1/3  in  the  24-stream 
model. 

The  DOM  using  Gaussian  quadrature  is  useful  when  dealing  with  RT  problems 
in  that  the  choice  of  directions  can  be  closely  matched  to  the  means  of  expressing 
the  scattering  phase  function  as  a  low  order  Legendre  expansion.  In  general,  the 
level  of  detail  in  a  discrete  ordinate  method  is  directly  related  to  the  level  of 
detail  desired  in  the  Legendre  expansion  of  the  phase  function.  That  is,  using 
the  discrete  ordinates  technique,  one  approximates  the  phase  function  (P(/i)) 
by  a  finite  sum  of  a  series  of  weighted  Legendre  polynomials  as 
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where  PpfJ,)  is  t'th  Legendre  polynomial. 
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and  /i  is  the  cosine  of  the  angle  between  incident  and  scattered  directions. 

For  a  given  discrete  ordinates  model  (Sjv),  the  limit  (N)  of  the  summation  P(/i) 
is  the  same  as  the  limit  of  the  stream  model  N  value  used.  This  ensures  the 
adequate  precision  of  angular  integrals.  However,  the  summation  P(/i)  is  seldom 
used  in  the  raw  form  seen  above.  Instead,  the  scaling  transformation  is  used  to 
modify  the  Legendre  polynomial  weighting  constants  Xi  to  remove  the  forward 
peak  of  the  phase  function.  The  appendix  describes  a  modified  version  of  the 
scaling  transformation  approach  that  optimizes  the  choice  of  Legendre  weighting 


22 


coefficients.  This  approach  also  results  in  changes  to  the  coefficients  a  and  ru 
used  in  the  scattering  computations  on  a  cell-by-cell  basis.  In  the  development 
that  follows,  this  step  has  been  assumed  completed  such  that  the  resulting  phase 
function  has  no  forward  peak. 

From  our  discrete  ordinates  approach,  then,  the  scattering  source  function 
for  propagation  direction  is 


Jmjr)  =  zu{r)  — ^  (21) 

m'  =  1 

2.4  Low-Density  DOMRT 

Armed  with  this  understanding  of  the  discrete  ordinates  approach  to  angular 
integration,  from  first  principles  we  proceed  to  generate  a  set  of  equations  for 
all  diffuse  streams  and  to  determine  the  net  radiation  flowing  in  any  direction 
at  any  position  within  the  scattering  volume. 

Let  s  be  a  distance  variable  in  the  direction  of  propagation  Then  Eq.  (15) 
can  be  written  as  a  function  of  distance  s  across  the  cell, 

dTp  (s) 

— +  (3)  =  0,  (22) 

which  has  the  explicit  solution, 

exp(-(7s),  (23) 

where  as  is  an  optical  depth  measure,  and  exp(— crs)  is  a  transmittance. 

Converting  back  to  the  original  variables,  we  find  a  solution  for  the  streaming 
radiance  over  a  single  cell, 

4,^ (■5)=  4,^(0)  exp(-(T5) +  Tp^  [1 -exp(-(7s)],  (24) 

where  Tp  =  Tp  (0)  =  Tp  (s)  is  the  sum  of  the  source  contributions  to 

'“m  '"m  '  ^  '"m  '  ^ 

radiance,  assumed  constant  within  the  cell  (c.f.,  Lewis  and  Miller  1984). 

2.4’ 1  Propagation  Equation  Terms 

Under  traditional  DOM,  the  transmission  factor  is  linearized,  whereby 

exp(  — crs)  sa  [s“^/(s“^  -|- cr )].  (25) 

Figure  1  compares  this  approximation  with  the  actual  negative  exponential 
function  via  a  plot  of  exp(— cr  s)/[s“^ /(s“^  -|-  cr)].  As  the  figure  illustrates, 
this  approximate  function  is  only  valid  for  very  small  optical  depths.  For 
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Figure  1.  Ratio  of  the  negative  exponential  function  exp(— crs)  to  the 
approximation  curve  +  '^  )]•  The  two  functions  are  only  roughly  equal 

around  an  optical  depth  of  zero.  (Optical  depth  is  r  =  as.) 


larger  optical  depths  the  approximate  transmittance  far  exceeds  the  actual 
transmittance  function. 

Introducing  Eq.  (25)  into  Eq.  (24)  produces  the  result, 


.-1 


+  CT 


(26) 


where  s  is  usually  replaced  by  the  distance  to  the  center  of  the  cell. 

In  standard  texts  describing  this  approach,  the  streaming  energy  is  discretized 
such  that  a  stream  directed  in  the  direction  at  position  (xp  j/j,  Zk)  is  assigned 
a  value  Stream  radiances  at  walls  are  identified  by  the  use  of  half 

indices.  Eor  example,  for  a  stream  flowing  into  the  first  octant,  there  would 
be  surface-based  stream  values  at  the  x  U  (-^(m,i,j-i/2,A;) )? 

-2^  {^(m,i,j,k-i/2))  input  faces.  On  each  of  these  input  faces,  the  modeled  net 
energy  crossing  the  face  is  proportional  to  fxAyAz  on  the  x  input  face,  rjA^Az 
on  the  y  input  face,  and  ^A^Ay  on  the  2;  input  face,  where  the  A’s  are  the  cell 
lengths  along  the  respective  subscripted  axes. 

Thus,  in  each  case,  the  net  energy  crossing  the  input  face  of  the  cell  is 
proportional  to  the  cosine  of  the  stream  with  the  input  face  multiplied  by  the 
surface  area  of  the  face  itself.  The  mean  source  radiance  at  the  cell  edge  can 
therefore  be  expressed  as  a  surface  average  over  the  three  input  surfaces: 


4.0) 


flAy  Azl(^m,i  — 1/2,  j,k)  T  ^^xAzI(^rn,i,j  —  l/2,k)  T  Ay  I(^zn,i,j,k  — 1/2) 


liAyAz  +  f]A,Az  +  CA,A, 


(27) 
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Dividing  numerator  and  denominator  by  AxAyAz/2  we  obtain, 

m,i  — 1/2,  j,k)  /  2T]I(^zn,i,j  —  1  /2,k)  /  2^I(^zn,i,j,k  —  l 

2/1  /  Ax  +  2ri  /  Ay  +  2^/ Az 

(28) 

(0)  thus  represents  a  weighted  average  of  the  projected  intensities  on  the 
three  input  faces  of  the  volume  of  interest.  This  mean  initial  value  for  the 
stream  is  the  first  unknown  on  the  RHS  of  Eq.  (26).  The  second  unknown  is  the 
choice  of  the  distance  parameter  s  measuring  the  distance  to  the  center  of  the 
cell.  In  the  standard  approach  s  is  normally  determined  as, 

1  _  2/i  2?7  2^ 

s  Aj,  Ay  Az' 

where  Ax/2ii  would  be  the  total  distance  travelled  by  a  ray  crossing  a  plane- 
parallel  slab  of  thickness  Ax/ 2  that  makes  a  cosine  of  /i  with  the  slab  normal. 

The  represents  the  volume  centered  limiting  path  radiance  based  on  the 

I{m,i,j,k)  values  obtained  from  a  previous  iteration  of  the  model,  the  direct 
illumination  for  that  cell,  and  the  cell  centered  temperature  T(j  j  yt)- 
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2.^.2  Low-Density  Limitations 

Though  variations  on  the  above  described  approach  have  been  investigated,  one 
essential  point  remains:  the  intensity  projected  for  the  opposite  side  of  the  cell  is 
typically  derived  via  extrapolation  of  its  value  at  the  cell  center.  One  particular 
choice  for  this  extrapolation  is  known  as  the  diamond  difference  technique.  To 
illustrate  this  technique,  let  /_  be  the  initial  value  of  the  stream  at  the  near  edge 
of  the  cell  (previously  (0)),  let  Jo  be  the  cell-centered  value  (/^  (■5))5  ^-nd 
let  /_|_  be  the  estimate  for  the  stream  at  the  far  side  of  the  cell  (/^  (2^)).  Under 
the  diamond  difference,  the  central  stream  value  is  simply  approximated  as  the 
mean  of  the  stream  values  at  the  opposite  sides  of  the  cell  (Jo  =  \{I-  +  -^-i-))- 
To  estimate  the  stream  value  at  the  far  side  of  the  cell,  this  equation  is  simply 
inverted  to  produce: 


/+  =  2/o  -  /. 


s  ^  I  —  tj  Tq 

2 - ^  -  1. 

s  ^  a 


^  ^  (2-^0  -  I-) 

-|-  u 


(30) 


Difficulties  arise  whenever  2L^  <  /_  because  the  numerator  has  the  potential 
to  become  negative  for  large  optical  depths.  For  example,  what  happens  when 
modeling  a  cloud  covering  more  than  one  cell?  Normally,  models  begin  by  setting 
all  radiances  initially  to  zero.  Hence,  on  the  first  pass,  =  0.  Figure  2  shows 
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Figure  2.  Negative  flux  consequences  of  the  diamond  difference  equation. 


the  behavior  of  in  this  case  as  a  increases.  This  indicates  that  for  cell  optical 
depths  greater  than  2,  the  flux  will  be  estimated  as  negative. 

This  problem  must  be  patched  in  order  to  hope  for  a  reasonable  solution.  One 
patch,  the  so-called  negative  flux  fixup,  yields  a  numerically  stable  algorithm 
whereby  all  negative  fluxes  are  simply  reset  to  zero.  However,  the  result  of 
using  this  fix  produces  difficulties  in  interpreting  the  appropriate  brightness  of 
regions  where  these  negative  values  persist  throughout  successive  iterations  of 
the  modeled  volume.  In  the  case  mentioned  above,  where  a  cloud  spans  more 
than  one  model  cell,  the  zero  values  produced  on  outer  cells  may  never  allow 
energy  to  propagate  into  the  volume  to  reach  an  inner  cell.  The  algorithm  then 
predicts  black  regions  that  are  physically  unreasonable  in  terms  of  typical  cloud 
scattering  properties.  Some  better  means  of  treating  optically  dense  clouds  is 
thus  sought  here. 

Other  authors  have  proposed  modifications  of  the  standard  discrete  ordinate 
method  (c.f.  Wakil  and  Sacadura  (1992)  and  Chai  et  al.  (1994)).  However,  these 
modifications  retain  the  usage  of  a  cell-centered  computation  of  the  limiting  path 
radiance  and  still  require  renormalization  to  produce  energy  balance.  In  the  next 
section,  a  theoretical  discussion  of  a  discrete  ordinates  radiative  transfer  solution 
is  presented  which  extends  the  traditional  volume-based  concepts  to  a  surface- 
based  method.  It  is  shown  that  this  reworking  of  the  discrete  ordinates  approach 
improves  the  performance  of  the  discrete  ordinates  results  for  dense  media  and 
can  support  a  system  which  automatically  conserves  energy. 
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2.5  Surface-Based  DOM 


In  this  section,  an  approach  is  developed  to  track  available  diffuse  energy 
sources  using  surface-based  input  and  output  stream  radiances  over  an  individual 
scattering  cell.  We  first  discuss  transmission  across  the  cell  between  input  and 
output  faces.  We  then  focus  on  evaluating  the  volume  averaged  unscattered 
radiation  due  to  energy  entering  the  cell  through  a  single  cell  wall  along  a 
single  streaming  direction.  Using  this  result,  and  assuming  the  diffuse  energy 
available  in  a  given  streaming  direction  is  constant  across  a  given  input  surface, 
one  can  determine  the  influence  of  that  stream  on  the  path  radiance  of  any 
stream  as  it  exits  the  volume  through  a  given  surface  element.  Also,  from  the 
standpoint  of  an  output  face,  the  average  transmission  through  the  volume  for 
a  given  stream  exiting  a  given  face  can  be  determined.  From  this  information, 
one  can  compute  the  average  contribution  of  volume-based  path  radiance  to  an 
output  stream  from  all  input  streams.  This  information  can  then  be  combined 
to  determine  the  net  fractional  energy  accounted  for  due  to  volume-based 
transmission,  absorption,  and  diffuse  scattering  processes  for  energy  entering 
a  given  face  via  a  given  stream.  It  is  argued  that  the  remaining  energy  can  be 
accounted  for  by  a  surface  scattering/absorption  procedure. 

Since  all  results  must  be  expressed  in  terms  of  output  surface  radiances,  this  is 
best  accomplished  by  averaging  Eq.  (24)  over  an  output  face, 

exp(-as)^  +  [1  -  exp(-(7  s)] ^  .  (31) 

Note  that  already  represents  a  volume  average;  what  remains  are 

evaluations  of  the  output  face  averaged  transmission  factor  and  the  expectation 
value  of  the  product  factor,  which  is  the  first  term  on  the  RHS  of  Eq.  (31).  In 
this  calculation,  we  are  aided  by  the  assumption  that  all  input  radiances  are 
constant  over  each  input  face.  This  first  term  thus  is  reexpressed  as  a  weighted 
sum  of  terms  as  explained  in  the  following  sections. 

^,5,1  Cell  Geometry 

In  describing  these  expectation  values,  we  begin  by  representing  the  directional 
vector  O  by  its  three  directional  cosines  (/i,?/,!^)  in  x,  y,  and  2;  directions, 
respectively.  Because  the  cells  have  a  cubical  geometry,  we  can  assign  a 
‘standard’  stream  chosen  to  flow  into  the  first  octant,  principally  in  the 
direction.  Hence,  ^  ^  >  t]  >  Q.  The  appearance  of  the  cell  with  respect 

to  the  incident  stream  is  illustrated  in  figure  3  as  it  appears  to  an  incident 
stream.  In  the  figure,  the  2;  axis  extends  into  the  page,  and  the  origin  is  labeled 
O  ,  near  the  lower  left  corner. 

Eor  any  stream,  there  will  always  be  three  input  faces  and  three  output  faces. 
We  identify  the  various  faces  of  the  volume  using  a  sign  (+/-)  and  letter  (X,  Y, 
Z)  convention.  In  the  standard  scenario,  because  1^,  y,  and  rj  are  all  greater  than 
zero,  the  input  faces  will  be  identified  as  the  -X,  -Y,  and  -Z  faces;  the  output 
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Figure  3.  Division  of  a  cell  for  transmission  purposes.  Face-on  view  as  seen  by 
energy  in  the  incident  stream. 


faces  will  be  the  -|-X,  -|-Y,  and  -\-Z  faces.  We  subdivide  each  of  the  3  input  faces 
into  a  set  of  rectangular  and  triangular  faceted  subareas  as  shown  in  figure  3. 
There  are  four  types  of  such  subareas,  labeled  with  letters  A,  B,  C,  and  D. 

To  better  understand  the  subdivision  of  the  cell  input  faces  into  these  different 
facets,  we  provide  two  additional  figures.  Figure  4  shows  the  cell  segmenting 
system  as  the  stream  enters  and  exits  the  cell  from  a  rotated  viewing  position. 
Figure  5  shows  the  same  stream  and  segmenting  system;  dashed  lines  have 
been  added  to  show  the  hidden  boundaries  of  the  segmented  portions  within 
the  volume.  In  addition,  letters  in  figure  5  denote  the  corners  of  the  volume 
segments. 

In  table  1,  the  descriptions  of  the  input  and  output  faces  of  each  volume  segment 
are  given  to  clarify  the  connection  between  figure  3  and  figures  4  and  5.  The 
designation  of  the  appropriate  input  and  output  faces  with  which  each  segment 
is  associated,  along  with  the  appropriate  corner  point  labels,  are  given  to  help 
understand  the  3D  nature  of  the  volume  and  its  segments. 

In  the  type  A  area,  the  optical  depth  is  constant  at  each  point  over  the 
rectangular  input  face.  In  figure  5,  this  corresponds  to  an  optical  depth 
experienced  over  the  path  from  input  point  a  to  output  point  n.  Geometrically, 
this  volume  is  a  parallelepiped.  For  type  B  and  C  input  regions,  the  geometry 
is  that  of  a  wedge:  the  optical  depth  is  constant  along  one  axis  and  linearly 
increasing  along  the  second.  Type  D  areas  represent  the  bases  of  trigonal 
pyramids,  with  a  maximum  optical  depth  along  the  edge  bordering  the  A  region, 
and  have  optical  depth  decreasing  linearly  down  to  zero  on  the  opposite  side. 


28 


Figure  4.  Subdivided  cell  revealing  input  and  output  surfaces  in  a  rotated  view 
for  the  same  stream  used  in  figure  3. 


P  q  r 


Incident  Stream  Direction 


Figure  5.  Subdivided  cell  showing  hidden  edges  bordering  different  volume 
segments  in  a  rotated  view. 
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Table  1.  Input  and  output  face  designations  for  surface  facets 


Region 

Input  Face 

Input  Surface 

Corners 

Output  Face 

Output  Surface 

Corners 

A 

-Z 

+z 

n,o,q,r 

Bi 

-X 

a,b,k,l 

+z 

k,l,n,o 

B2 

-Z 

d,e,g,h 

+x 

g,h,q,r 

Cl 

-Y 

a,d,m,q 

+z 

m,n,p,q 

C2 

-Z 

b,c,e,f 

+Y 

c,f,o,r 

Di 

-Y 

a,j,m 

+z 

j,m,n 

D2 

-X 

a,j,k 

+z 

j,k,n 

Ds 

-Z 

e,h,i 

+x 

h,i,r 

Di 

-Z 

e,f,i 

+Y 

f,i,r 

D5 

-Y 

d,g,p 

+x 

g,p,q 

De 

-X 

b,c,l 

+Y 

c,l,o 

2.5.2  Facet  Transmittance s 

From  the  standpoint  of  the  discrete  ordinates  method,  all  the  energy  in  a  given 
stream  is  considered  to  be  traveling  in  the  same  direction  (O^),  even  though  it  is 
applied  over  a  specific  solid  angle  during  integration.  Second,  the  geometry  of  all 
cells  is  assumed  cubic.  Thus,  every  cell  in  the  scattering  volume  is  characterized 
by  a  single  cell- width  parameter  A.  Energy  entering  a  cube  through  the  area 
labeled  A  experiences  a  total  optical  depth,  r,  equal  to  crA/i^,  as  it  crosses  the 
volume.  So  the  area-averaged  mean  transmittance  for  flux  entering  the  volume 
through  area  A  is 

Ta  =  exp(-r).  (32) 

Area  types  B,  C,  and  D  have  more  complicated  forms  for  the  transmittance 
equations  due  to  their  modified  geometries. 

For  type  B  and  C  areas,  the  mean  transmittance  can  be  computed  by  taking 
the  integral  of  the  point-to-point  transmittance  at  each  point  on  the  surface  and 
dividing  by  the  total  surface  area  to  produce: 

1 

-  -  f  e~'^  -  1 

Tb  =  Tc  =  /  e'x\){—Tu)  du  =  - .  (33) 

0 

Notice  that  this  result  depends  on  r  =  crA/i^,  but  not  on  the  exact  length  and 
width  of  the  region,  since  it  is  surface  averaged. 
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For  type  D  areas,  the  mean  transmittance  is  again  determined  by  integration: 

1  U 


Td  =  2  du  dv  exp(—Tv)  = 


-  (1  -  ^) 
t^I2 


(34) 


0  0 

Here,  division  by  the  area  normalization  factor  (1/2)  results  in  the  factor  2 
appearing  in  front  of  the  integral. 


For  the  sake  of  brevity,  these  results  can  be  described  using  an  apparently  new 
class  of  related  functions: 


G*(r)  = 


oo 


(— r)*/z! 


(-T)' 


(35) 


Thus,  we  can  represent  the  average  transmittances  for  cell  wall  regions  A,  B,  C, 
and  D  as  Ta  =  Tb  =  Tc  =  Gi(r),  and  To  =  G2{t).  Comparisons  of 

the  G  functions  are  given  in  figures  6  and  7.  The  first  of  these  illustrates  the 
basic  behaviors  of  these  functions.  (We  extend  these  results  up  to  ^3(7")  because 
this  functional  form  is  relevant  in  further  calculations.)  Figure  7  illustrates  the 
normalized  behaviors  of  this  class  of  functions.  Note  that  successive  functions 
have  limiting  behaviors  Giir)  sa  1  —  r/(z  +  1)  for  small  r.  But,  even  after 
normalizing  these  functions  according  to  their  limiting  behavior  around  zero, 
higher  order  G  functions  are  always  greater  than  lower  order  functions.  This 
illustrates  the  fact  that  energy  passing  through  a  cell  between  adjacent  walls 
experiences  far  less  attenuation  than  is  predicted  by  the  negative  exponential 
function,  and  partially  explains  why  the  classical  FE  DOMRT  methods  fail. 


Optical  Depth  (r) 


Figure  6.  Comparison  of  transmittance  functions  Giir)  for  i  values  ranging  from 
0  to  3. 
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Figure  7.  Comparison  of  transmittance  functions  6i(r)  normalized  according 
to  their  limiting  behavior  around  r  =  0. 


2.5.3  Face-to-Face  Transmittance s 

Using  these  results,  one  may  describe  the  mean  transmittance  of  energy  flowing 
from  one  input  face  to  a  given  output  face  through  weighted  averages  of  the 
transmittances  of  its  different  facets. 

To  produce  these  transmission  functions,  the  parameters  a  =  <  1,  and 

P  —  I7/CI  ^  7  ^re  introduced.  We  also  introduce  area  measurements  that  are 
determined  via  projections  onto  a  plane  containing  the  -Z  face  (figure  3).  The 
projected  areas  subtended  by  the  various  surface  portions  with  respect  to  this 
plane  are  given  in  table  2. 

Table  2.  Region  area  values  for  subcomponents  of  input  surfaces 


Region  Type 

Region  Area 

A 

B 

a(l  -  /3)A2 

C 

(1  -  a)/3A2 

D 

a/3AV2 

One  can  then  determine  the  average  transmittance  from  a  given  input  wall  to 
a  given  output  wall  by  forming  a  weighted  average  of  the  transmittances  of 
the  various  facets  of  each  input  wall  connecting  to  each  output  wall.  These 
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calculations  form  a  matrix  with  components  detailed  as  follows:  The  -Z  and  +Z 
walls  only  connect  through  the  type  A  area.  Hence, 


T-z,+z  =  Gq{t  ). 


(3C) 


The  X  and  Z  walls  connect  via  a  D-type  and  a  B-type  volume  segment.  Let  As 
and  Ad  be  the  projected  areas  of  these  walls  respectively.  Thus, 


T_ 


z,+x 


=  T_ 


x,+z  = 


Ab  Tb  +  Ad  Td 
Ab  +  Ad 


(37) 


which  results  in. 


T_ 


z,+x 


=  T_ 


x,+z  = 


2Gi{t)  —  /3[2Gi(r)  —  G2{t)\ 
2-/3  ' 


(38) 


where  the  a  factor  cancels  out  of  the  numerator  and  denominator.  A  similar 
development  can  be  performed  for  the  connectivity  between  Z  and  Y-type 
surfaces, 

rp  rp  2Gi{t)  -  a[2Gi{T)  -  G2{t)] 

J--Z,  +  Y  =  J--Y,  +  Z  =  - V - •  [oy) 

I  —  a 

The  X  and  Y  faces  only  connect  through  single  facets  (D5  and  De)  and  thus  are 
characterized  by, 

T-x,+y  =  T-y,+x  =  G2{t  ).  (40) 

Lastly,  neither  the  X  nor  Y  faces  connect  directly  between  their  input  and  output 
faces: 

T-y,+y  =  T-x,+x  =  0.  (41) 

Subscripts  in  these  equations  indicate  the  face  connections.  The  resulting  matrix 
is  diagonally  symmetric. 

In  addition  to  the  foregoing  transmittance  factors  connecting  input  and 
output  surfaces  for  the  standard  scenario,  we  may  also  determine  the  mean 
transmittance  factors  to  each  output  face  (+X,  +Y,  and  +Z)  to  all  connected 
input  faces.  These  calculations  are  accomplished  by  constructing  a  weighted 
sum  of  the  previous  transmittance  factors,  where  the  weights  are  determined 
according  to  the  fraction  of  the  output  face  connecting  to  each  input  face.  For 
an  X-type  output  face,  the  energy  is  passing  through  2  D-type  facets  and  1 
B-type  facet.  Thus, 

Ab  Tb  +  2  Ad  Td 

Tx  =  °  d  „  .  (42) 


Ab  +  2  Ad 


which  results  in. 


fx  =  Gi(r)-/3[Gi(r)-G2(r)]. 
A  similar  result  is  obtained  for  Y-type  faces, 

Ty  =  G,{t)  -  a[G,{T)  -  G2{t)]. 


(43) 

(44) 
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For  Z-type  faces,  the  summation  averages  the  results  over  an  A-,  a  B-,  a  C-,  and 
two  D-type  facets,  leading  to  the  result, 

fz  =  ^8Gq{t)  +  {a8  +  7/3)Gi(r)  +  a/3G2('r),  (45) 


where  7  =  1  —  a  and  8  =  1  —  fi. 

Due  to  the  symmetries  of  the  problem,  Eqs.  (43)  through  (45)  also  represent  the 
average  transmittance  through  the  volume  from  a  given  input  face  of  the  same 
type.  That  is,  Tz  is  both  the  mean  transmittance  of  energy  entering  the  volume 
through  a  Z-type  input  face  and  the  average  transmittance  of  energy  exiting  the 
volume  through  a  Z-type  output  face. 

These  computations  allow  us  to  determine  both  the  total  amount  of  flux  lost  due 
to  scattering  and  absorption  from  the  stream  as  it  enters  a  given  face,  and  also 
allow  for  an  emission  term  to  be  evaluated  on  a  similarly  labeled  output  face. 
For  example,  if  the  average  transmission  through  the  volume  to  a  given  output 
face  is  labeled  Tq,  then  (1  —  Tq)  represents  a  source  factor  in  the  emittance 
portion  of  the  radiative  transfer  solution  for  that  output  face. 

Therefore,  let  us  return  to  Eq.  (31),  where  as  we  already  noted,  the  factor 
already  represents  a  volume  average  result.  Thus,  this  factor  may  be  removed 
from  the  expectation  operator,  resulting  in, 

=  (4,^(0)  exp(-(7s)^  “  (exp(-(7s))].  (46) 

The  quantity  within  expectation  operator  in  the  second  term  on  the  RHS  of  the 
equation  is  just  the  Tq  terms  we  have  just  derived. 

Eurther,  with  our  current  understanding  of  the  transmission  factors  between 
different  cell  faces,  we  can  evaluate  the  first  term  on  the  RHS  directly.  Recall 
that  one  of  our  assumptions  was  that  the  input  radiances  were  assumed  constant 
over  each  input  face.  We  therefore  introduce  the  terms  /x,  /y,  and  Iz  for  our 
standard  scenario  X,  Y,  and  Z  input  face  radiances,  respectively.  The  first  term 
on  the  right,  under  our  constant  input  assumption,  becomes  a  weighted  sum 
of  the  area  of  a  given  facet  times  the  input  radiance  applicable  for  that  facet 
times  the  average  transmission  factor  for  that  facet.  Eor  example,  for  an  X-type 
output  face  we  have. 


(0)  exp(  — cr  s)'^ 


Ao  To  Iy  +  {  Ab  Tb  +  Ao  To)  Iz 


X  Ao  +  Ab  +  Ao 

2  V  2 


T_ 


Z,YX 


Iz, 


(47) 


where  the  quantities  /3/2  and  (1  — /3/2)  are  fractional  area  weighting  factors  that 
sum  to  unity  and  represent  the  contributions  of  the  input  faces  to  the  output 
from  a  given  output  face.  The  weights  for  each  output  face  are  given  as  a 
function  of  the  appropriate  standard  scenario  input  faces  in  table  3.  Of  course, 
transformations  must  be  made  to  translate  these  results  into  terms  useable  for 
non-standard  scenarios. 
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Table  3.  Input  face  weight  values  for  the  standard  scenario 


Output 
Face  Type 

X  Input 
Weight 

Y  Input 
Weight 

Z  Input 
Weight 

X 

0 

/3/2 

(l-,8/2) 

Y 

a/2 

0 

(1-0/2) 

Z 

a(l-f3/2) 

[3(1 -a/2) 

1 

1 - 1 

1 

1 - 1 

2.5.4  Mean  Illumination  Factors 

We  have  thus  far  developed  a  mechanism  for  evaluating  all  the  terms  in  the 
average  equation  of  transfer  between  input  and  output  faces  except  for  the 
volume- averaged  limiting  path  radiance,  .  The  limiting  path  radiance  {L) 
is  the  sum  of  blackbody  (-B),  scattering  source  (J),  and  direct  contributions.  In 
considering  the  computations  needed  to  evaluate  the  radiative  transfer  equation, 
we  first  acknowledge  that  J,  as  given  in  Eq.  (10),  is  a  positionally  varying 
quantity.  But,  I  is  not  known  except  as  wall-average  values.  We  are  thus 
limited  according  to  the  granularity  of  the  volume  resolution  used.  Further, 
since  the  phase  function  reflects  the  results  of  single  scattering,  rather  than 
posit  a  complicated  model  for  /’s  behavior  within  the  volume,  to  first  order 
it  is  postulated  that  the  scattering  due  to  a  given  input  stream  is  dependent 
on  only  the  mean  unscattered  illumination  provided  by  that  stream  and  the 
medium  within  a  given  cell.  The  determination  of  this  mean  value  requires  a 
set  of  integrals  similar  to  those  used  to  determine  the  mean  transmittances. 
These  quantities  are  produced  by  integrating  over  the  scattering  volume  swept 
by  energy  passing  through  geometries  identical  to  those  considered  during  the 
transmittance  discussion. 

For  each  type  region,  the  unscattered  illumination  present  at  a  certain  depth 
within  the  material  is  equal  to  exp(— where  r  is  the  maximum  optical 
depth  across  the  cell  for  that  material  and  for  that  stream,  and  where  w  is  the 
fractional  distance  of  that  point  into  the  material  compared  with  the  maximum 
optical  depth  presented  by  that  material  for  that  stream. 

For  region  A,  the  mean  illumination  can  be  related  to  an  integral  over  a  cubic 
region  where  the  optical  depth  is  measured  into  the  material  along  a  single 
axis.  A  unit  volume  is  used  because  the  optical  depth  along  the  w  axis  is 
measured  such  that  at  distance  1  the  maximum  optical  depth  is  reached.  Along 
the  other  axes,  the  averaging  process  eliminates  any  dependence  on  actual  region 
size.  Additionally,  this  geometric  integral  is  possible  due  to  the  assumption  of 
uniform  initial  radiance  of  the  incident  stream  across  the  input  face.  These 
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considerations  lead  to  the  formula, 


A  = 


111 

J  du  J  dv  J  exp(—Tw)  dw 

0  0  0 _ 

111 
J  du  J  dv  J  dw 
0  0  0 


e-^  -  1 


— r 


=  G.(r). 


(48) 


For  B  and  C  type  areas,  the  following  integral  can  be  used.  In  this  integration, 
the  u  axis  denotes  the  direction  of  decreasing  wedge  thickness. 


B  =  T  c  = 


1  1  l-u 

J  du  J  dv  J  exp(—Tw)dw 

0  0  0 _ 

1  1  l-u 

J  du  J  dv  J  dw 
0  0  0 


e  ^-(1-r) 


72 


=  G2(r),  (49) 


For  type  D  areas,  the  integral  is  performed  over  a  pyramidal  shaped  volume. 
The  thickest  portion  of  this  region  is  at  u  =  0,  u  =  0  in  the  following  integral. 


D  = 


1  U  1—u 

J  du  J  dv  J  exp(—Tw)  dw 

0  0  0 _ 

1  U  1—U 

J  du  J  dv  J  dw 
0  0  0 


-  (1  -  r  +  r^/2) 
—t^/Q 


=  Gs(t),  (50) 


These  ,  functions  are  unitless  and  relate  to  the  relative  mean  illumination 
(radiance)  provided  by  flux  entering  the  volume  segment  of  interest.  To  produce 
overall  average  illumination  in  the  cell  due  to  flux  flowing  across  a  given  input 
face,  the  ‘standard’  model  is  again  referred  to.  We  denote  the  relevant  volumetric 
mean  illuminations  due  to  incident  radiances  /x,  /y,  and  Iz  (passing  through 
the  standard  input  faces  X,  Y,  and  Z  of  the  volume,  respectively)  by  the 
quantities  /x,  /y,  and  Iz-  These  means  are  computed  based  on  weighted 
averages  of  the  mean  illumination  of  the  volume  type  illuminated  [Eqs.  (48) 
through  (50)],  multiplied  by  the  fraction  of  the  total  volume  contained  within 
each  volume  segment.  The  segment  volumes  associated  with  each  of  the  region 
types  are  given  in  table  4. 

Table  4.  Segment  volumes  for  subcomponents  of  input  region  types 


Region  Type 

Region  Volume 

A 

1 

1 — 1 

1 

1 — 1 

B 

a{l  -  /3)AV2 

C 

(1  -  a)/3AV2 

D 

q;/3A^/3 
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For  the  resulting  mean  illuminations,  one  then  obtains  (after  dividing  by  the 
total  volume  A^): 


T7=.xIx=abhA+,3^}/v,  (51) 

fr  =  ,v-Iy=;3{T^  +  ah^}ly,  (52) 

and 

Iz  =  ■,  z^z  =  |7(^Gi(r)  +  [/37  +  a6]  (53) 

The  total  mean  illumination  within  the  scattering  volume  for  a  given  stream  is 
then  just  the  sum  of  the  mean  illuminations  over  the  appropriate  three  input 
faces  for  that  stream.  Due  to  the  symmetrical  structure  of  the  streams  in 
different  octants,  it  is  possible  to  generate  ,  x,  ,  Y,  and  ,  z  values  according  to 
stream  type  and  make  appropriate  transformations  to  convert  the  input  stream 
information  from  the  proper  input  faces  into  the  mean  illumination  statistics  for 
the  cell. 

Based  on  the  mean  illumination  for  a  given  stream,  we  may  then  determine 
the  amount  of  flux  subsequently  scattered  into  each  stream  leaving  the  volume 
from  each  face,  and  thus  assess  the  total  flux  leaving  the  cell  due  to  volumetric 
scattering. 

2.5.5  8-Stream  Example 

To  illustrate  the  development  to  this  point,  consider  an  8-stream  DOM  model. 
In  this  model,  |/i|  =  \ri\  =  |i^|  =  1/a/3-  Hence,  a  =  jS  =  1.  Because  of  model 
symmetries,  all  eight  streams  can  be  characterized  by  transformations  on  the 
one  unique  standard  stream  direction  in  the  first  octant. 

As  each  input  stream  enters  through  a  particular  input  face,  a  portion  of  the 
entering  flux  will  transmit  directly  to  two  output  faces,  to  each  of  which  half  the 
streaming  flux  is  directed.  The  transmittance  to  each  output  face  is  G2(t),  where 
r  =  \/3(7A.  On  each  output  face,  Q,  four  streams  are  directed  outward  through 
it.  Each  of  these  output  streams  has  a  mean  emittance  factor  of  (1  —  Tq)  = 
[1  —  G2{t)].  Thus,  the  single  scattered  flux  emitted  by  a  particular  stream  on  a 
given  face  is  proportional  to  [1  —  G2{t)]. 

For  each  stream  m,  the  volume  averaged  radiance  is 

T  \  ^  ^ )  f  K  A  \ 

Im  =  Gsir) - - - .  (54) 

Each  input  stream  scatters  flux  into  all  eight  output  streams.  Each  stream 
makes  a  cosine  of  1  with  itself,  -1  with  the  stream  in  the  opposite  octant,  1/3 
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with  three  streams,  and  -1/3  with  the  three  remaining  streams.  Let  these  cosines 
be  represented  by  the  variable  (cosine  of  the  scattering  angle),  and  represent 
the  modeled  phase  function  value  for  this  scattering  angle  by  For  the 

moment  ignoring  emissive  effects,  the  influence  of  input  stream  m  on  a  specific 
outward  flowing  stream  s  is  computed  using  the  volume  averaged  radiance  and 
the  phase  function  value.  Now,  if  we  designate  the  scattering  source  in  this 
direction  as  Jms-,  then  from  Eq.  (21)  we  have, 

47r?x>  ~  — 

Jms  —  Z  P{pms  )  J m  ■  (55) 

o 

Here,  the  angular  weight  factor  w  equals  1  for  all  8-stream  cases.  This  source 
term  can  then  be  used  to  determine  input  stream  m’s  effect  on  all  scattered  flux 
exiting  the  three  output  faces  corresponding  to  stream  s.  Let  E^s  be  the  net 
flux  exiting  the  volume  due  to  Jms- 

1  A-7nn 

=  3^^[1  -  G2{t)]  (56) 

v3  o 

where  the  3  represents  the  contributions  on  the  three  output  faces,  1/a/3  is  each 
output  stream’s  cosine  with  each  output  face  normal,  47^^/;/8  appears  due  to  the 
integration  over  solid  angle  on  the  output  face,  and  [1  —  G2{t)]  is  the  emittance 
factor.  The  units  of  E^s  are  W / fiva.  Summing  over  all  exiting  streams,  the 
total  flux  streaming  out  of  the  volume  due  to  scattering  source  J^s  is: 

Es  =  Y,E^,  =  ^^[1-  G2{r)]  ^Gs{t)7;;Y.Y 

Due  to  the  discrete  ordinates  method  chosen,  the  sum  over  scattering  angles 
can  be  divided  into  individual  sums  over  each  Legendre  component.  These 
sums  cancel  for  all  but  the  isotropic  term,  which  sums  to  unity.  Thus  through 
cancellation  we  obtain: 

1  A-tt 

Es  =  wA^——[1-  G2{r)]  Gs{t)  {Ix  +  Iy  +  Iz)-  (58) 

v3  o 

Similar  computations  can  be  made  for  the  total  transmitted  {Et)  and  total 
incident  fluxes  (Ej)  due  to  this  stream  (again  taking  account  for  the  finite  solid 
angles  associated  with  the  streams  as  they  enter  and  exit  the  volume): 

1  Ayr 

Et  =  A^  —  G2{t )  {Ix  +  Iy  +  Iz),  (59) 

and 

1  47r 

El  =  A^  —  {Ix  +  Iy  +  Iz)-  (60) 
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To  evaluate  the  fractional  flux  accounted  for  via  transmittance  and  volumetric 
scattering  and  absorption,  Es  is  divided  by  the  single  scattering  albedo  to 
account  for  absorbed  and  scattered  radiation.  This  result  is  added  to  Et,  and 
the  sum  is  divided  by  the  incident  flux.  The  result  is: 

A(r)  =  =  G,(t)  +  G3(t)  -  G3(r)G3(T) 

=  1  -  [1  -  G2{t)]  [1  -  G3(r)]  . 

This  function  is  plotted  in  figures  8  and  9.  These  figures  show  the  behavior  of 
A(r)  as  a  function  of  increasing  cube  axial  optical  depth  (crA).  From  the  plots, 
it  is  clear  that  the  volume  based  equations  are  reasonably  accurate  for  r  ^  1, 
since  in  the  limit  of  r  =  0  the  curve  has  slope  zero  and  equals  unity,  accounting 
for  all  the  flux.  For  large  r  the  approximation  is  progressively  more  lossy. 


Figure  8.  Fractional  flux  accounted  for  in  8-stream  model  using  volumetric 
effects. 


The  main  point  to  keep  in  mind,  though,  is  not  that  these  calculations  indicate 
that  there  are  problems  for  high  optical  depth  -  indeed  there  are  problems  -  but 
these  are  problems  that  are  common  to  every  strictly  volume-based  equation 
set.  The  conclusion  is  that  volume-based  equations  alone  cannot  account  for  all 
the  scattered  and  transmitted  energy. 

Cases  of  higher  order  quadrature  sets  exhibit  similar  behavior  to  the  8-stream 
curve  shown  here,  since  the  limitation  is  not  in  the  number  of  streams  but  due 
to  the  single  scattering  approach  taken  to  model  the  scattering  phenomenon. 
Nevertheless,  these  results  can  be  shown  to  preserve  more  accuracy  than  the  low- 
density  approach  since  the  transmittance  factors  are  at  least  being  computed 


39 


1.0 


-  8-Stream  Model 


Figure  9.  Fractional  flux  accounted  for  in  8-stream  model  using  volumetric 
effects.  Plot  of  high  optical  depths. 


exactly,  and  we  are  making  our  single  scattering  calculation  based  on  a  true 
volume  average  of  the  radiance  for  each  stream.  The  low-density  equations  fail 
by  approximately  1  optical  depth.  These  equations  perform  well  up  to  and 
exceeding  1  optical  depth,  but  they  do  not  fair  well  beyond  about  4  optical 
depths. 

Hence,  the  most  critical  factor  limiting  the  use  of  the  pure  volumetric  approach 
obtained  thus  far  is  due  to  the  nature  of  the  atmosphere  itself:  at  visible 
wavelengths  natural  cloud  optical  depths  are  often  hundreds  of  optical  depths 
thick.  For  a  3D-RT  simulation  of  an  atmospheric  volume  containing  even  a 
simple  cumulus  cloud,  the  method  described  to  this  point  would  still  require 
millions  of  cells  to  adequately  maintain  flux  balance  via  purely  volumetric  effects. 

2.5.6  Surfacelike  Interactions 

To  augment  the  technique  developed  thus  far  enabling  it  to  treat  large  optical 
depth  conditions,  it  is  postulated  that  any  missing  energy  can  be  accounted 
for  by  using  surface-like  scattering  and  absorption  events.  This  hypothesis  can 
be  justified  based  on  the  following  two  empirical  arguments.  First,  consider  a 
stream  of  energy  directed  perpendicular  to  a  face  of  entry  into  a  volume  element. 
At  high  optical  depths,  negligible  energy  propagates  unscattered  through  the 
volume,  and  the  mean  unscattered  radiance  rapidly  approaches  zero  as  one 
moves  further  into  the  volume.  Physically,  this  means  that  virtually  all  the 
energy  is  interacting  with  the  media  within  a  short  distance  of  the  input  face. 
This  scenario  would  produce  near-zero  transmission  and  illumination  factors 
(resulting  in  a  near-zero  volumetric  efficiency  factor  A),  and  virtually  all  the 
energy  would  be  scattered  or  absorbed  at  or  near  the  surface.  In  a  second 


40 


thought  experiment,  consider  a  stream  that  enters  the  volume  at  an  incidence 
angle  which  is  nearly  tangent  to  the  input  face.  In  this  case,  we  do  not  care  how 
thick  the  media  is,  but  only  consider  the  cell  averaged  illumination  produced. 
Because  the  cosine  of  the  incident  stream  is  so  small,  the  fractional  volume 
illuminated  by  the  stream  is  very  small,  and  thus  whatever  scattering  does  occur 
must  occur  within  a  short  distance  of  the  entry  face.  In  either  case,  then,  a 
volumetric  scattering  model  is  inappropriate.  Thus,  it  is  argued  that  any  energy 
not  seen  when  accounting  for  transmittance  and  volumetric  diffuse  scattering 
and  absorption  can  be  modeled  with  a  surfacelike  scattering/absorption  process. 

To  characterize  the  amount  of  energy  to  be  recovered  due  to  surface  effects 
in  the  8-stream  example,  we  would  need  to  allow  fractional  energy  in  the 
amount  [1  —  G2{t)]  [1  —  ^3(7")]  to  interact  at  the  surface  of  the  medium  for 
each  incident  stream.  Let  us  characterize  this  fractional  surface  interaction  with 
the  parameter  S(t)  =  [1  —  A(r)].  For  stream  models  beyond  S2  the  form  of  these 
S  expressions  will  also  depend  on  the  stream  index  and  the  particular  input  face 
being  considered.  However,  for  the  8-stream  model  we  have  a  simple  expression 
that  is  the  same  for  all  streams  and  all  faces.  For  a  more  complicated  scattering 
model,  we  turn  to  the  next  most  complicated  case,  the  24-stream  S4  quadrature. 

2.5.7  24- Stream  Example 

For  the  S4  case  using  the  LSE  quadrature  (Fiveland  (1991);  Lewis  and 
Miller  (1984)),  each  stream  is  composed  of  two  components  with  magnitudes 
of  a  =  0.3500212,  each,  and  one  component  with  magnitude  b  =  0.8688903. 
Thus,  a  =  a/6  =  0.4028370  =  /3.  This  results  in  a  slightly  more  complex 
equation  set  than  produced  for  the  S2  set.  Let  us  define  7  =  1  —  a  =  0.5971630. 
Then,  (suppressing  G  function  dependencies  on  r). 


I  m  — 


G2  2  ^3 

+  “  T 


fIx+lY)  + 


Iz.  (62) 


The  scattering  source  term  (Jms)  for  this  case  is  the  same  as  given  in  Eq.  (55), 
except  that  w  =  1/3  for  each  stream  under  S4.  It  then  remains  to  determine 
the  energy  exiting  all  three  output  faces  for  this  stream.  Notice  that  two  of  the 
output  faces  have  stream  cosines  of  a  and  a  mean  transmittance  value  that  is 
different  from  that  for  the  third  output  face.  These  considerations  lead  to  the 
expression. 

Eras  =  b^  L{t)w  ^  (63) 

6  6 

where 

L(t)  =  1  -|-  2a  —  7^Go(t)  —  4a7Gi(r)  —  3a‘^G2(T).  (64) 

Due  to  the  symmetries  of  the  S4  case,  the  scattered  energy  as  a  function  of 
direction  only  depends  on  I rm  since  the  remaining  factors  are  only 

dependent  on  r.  Thus,  again,  E^s  is  summed  over  all  scattering  directions, 
and  as  before,  only  the  zero-th  order  Legendre  polynomial  has  a  nonzero 
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contribution,  such  that  the  total  volumetric  diffuse  scattering  produces  an 
output  energy  of, 

Es  =  —  L{t)  w  I  =  A^b  —  \L  (r)  ruT]  ,  (65) 

6  6  ^ 

which  depends  on  only  the  incident  energy  and  the  optical  depth  of  the  medium. 
Similarly,  the  incident  energy  and  the  transmitted  energies  only  depend  on  the 
optical  depth,  as, 

Ej  =  A^b^  [a  (Ix  +Iy)  +  Iz],  (66) 

6 

Et  =  AH^  [Tair)  a  {Ix  +  ly)  +  nir)  Iz]  ,  (67) 

6 

where 

Tair)  =  jGi{t)  +  aG2{T);  (68) 

Tbir)  =  j^Go{t)  +  2a7Gi(r)  +  a^G2{T).  (69) 

Combining  terms,  the  net  energy  unaccounted  for  via  transmission  and 
volumetric  scattering  can  be  expressed  as, 

Ei-Et-  Es/w  =  AH^[a  {Ix  +  /y)  [1  -  Aa(r)]  +  /y  [1  -  A6(r)]]  ,  (70) 

6 

Aa{T)  =  L{T)  +  +Ta{r),  (71) 

Ab{T)  =  L{t)  7^Gi(r)  +  2(^7^^^  +  +Tb{T).  (72) 

The  terms  Aa(r)  and  Ab{T)  represent  the  fractional  efficiencies  of  the  scattering 
model,  as  plotted  in  figure  10.  Standard  input  X-  and  K-type  faces  (those  where 
the  cosine  of  incidence  is  a)  are  handled  using  5*0  =  1  —  Aa(r),  and  Z-type  input 
faces  (those  with  a  cosine  of  incidence  of  b)  are  treated  using  S';,  =  1  —  A6(r), 
where  these  S  factors  are  used  as  the  fractional  amount  of  energy  arriving  at 
a  surface  that  must  be  scattered/absorbed  in  a  surfacelike  event  to  maintain 
energy  conservation. 

The  different  scattering  efficiencies  arise  because  the  cosine  of  incidence  results 
in  different  mean  path  lengths  through  the  media  depending  on  the  surface 
of  entry.  For  the  24-stream  case,  one  of  two  angles  of  incidence  are  possible. 
For  the  ‘6’  case,  there  is  a  greater  average  illumination  within  the  volume,  but 
this  is  offset  by  the  lower  average  transmittance.  The  net  result  is  that  Ab{T) 
has  a  uniformly  lower  efficiency  than  the  result  obtained  for  the  S2  order  case. 
This  result  is  obtained  even  after  accounting  for  the  larger  r  in  the  S2  case 
{t2  =  V^aA  versus  T4  =  crA/0. 8688903). 

Interestingly,  the  Aa(r)  curve  resembles  the  behavior  of  the  single  S2  case  curve. 
Perhaps  because  the  A^  curve  represents  a  smaller  cosine  with  the  entry  surface, 
while  the  transmittance  is  increased,  the  path  radiance  is  reduced,  and  the 
overall  performance  roughly  tracks  that  of  the  S2  order  case.  The  net  effect  is 
a  similar  behavior  (in  terms  of  overall  efficiency)  of  the  24-stream  case  to  the 
8-stream  case. 
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Figure  10.  Fractional  efficiencies  of  the  24-stream  model  using  volumetric  effects. 
The  ‘A’  curve  represents  the  function,  while  the  ‘B’  curve  is  the  A;,  function. 


2.5.8  ^8-Stream  Case  and  Beyond 

Beyond  the  24-stream  case,  there  are  increasing  numbers  of  equations  required 
to  account  for  the  efficiencies  of  each  scattering  stream.  For  the  48-stream 
approximation  there  will  be  four  independent  accounting  expressions:  two  for 
each  of  the  two  unique  weighting  structures.  In  the  80-stream  case  there  are 
three  unique  weighting  structures.  The  first  of  these  has  two  unique  curves; 
the  second  has  three  unique  curves  since  each  component  is  a  different  weight; 
and  in  the  third  case  there  is  one  unique  curve,  since  all  components  are  equal. 
There  are  thus  a  total  of  six  formulas  to  evaluate.  Obviously,  any  economies  to 
be  gained  through  symmetry  would  be  rather  limited  compared  to  the  added 
complexity  of  the  handling  routines  once  a  high  order  Legendre  expansion  model 
is  chosen. 

2.5.9  Accounting  for  the  Missing  Flux 

Given  that  the  means  of  accounting  for  missing  flux  is  through  surfacelike 
interactions,  how  should  one  approach  the  problem  of  correcting  the  radiance 
values  to  restore  the  missing  flux?  As  a  sample  case,  let  us  consider  radiance 
arriving  at  an  X-type  interface  in  the  24-stream  case  (cosine  of  incidence  with 
the  surface  normal  is  ji^n  =  o). 

We  have  already  decided  that  the  amount  of  flux  to  be  scattered  at  such 
a  surface  {Esx,m)  should  be  proportional  to  5'a(r)  and  the  incident  flux. 
We  also  know  that  the  scattering  process  itself  is  proportional  to  the  single 
scattering  phase  function  since  we  are  working  in  a  single  scattering 

paradigm.  And  recall  that  and  are  the  incident  and  scattered  radiance 
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propagation  vectors,  respectively.  Also  we  assume  that  Esx,m  can  be  related 
to  some  fraction,  of  the  overall  radiance  falling  on  that  surface.  These 

considerations  lead  to  the  equation: 

EsX,m  =  - l/^ml  dx,m  (73) 


where  Ix,m  is  the  mean  incident  radiance  at  that  face  for  stream  m. 

Since  we  assume  the  surface  scattered  radiance  {Iss,s,m)  at  this  interface  equals 


IsS,s,m  —  ^  dx,m  Ix,m  (4vr?X>j72 /8)  F(l2rns)i 
then  the  total  flux  scattered  at  the  surface  due  to  stream  m  will  be. 


E 


_  a2  4:7rWa  I  I  r  _  Esx,m  I 

SS,m  —  4a  y  ^  \l-ls  \  -I- SS,s,m  —  ^  m 

O  l/^m  I 


*1  r. 


(74) 


(76) 


One  can  easily  verify  that  this  result  actually  does  disperse  the  correct  amount 
of  flux  in  the  8-stream  conservative  isotropic  scattering  case  (where  P{jj.ms)  = 
(dvr)”^ )  case.  Here  we  have  w^n  =  Wg  =  1,  \fim  \  =  If^s  |  =  /i,  'n7  =  l,s  =  1...8. 
The  result  of  the  summation  in  Eq.  (75)  is  then  just  /i,  and  the  net  equation 
shows  that  Ess,m  =  Esx,m- 

In  general,  the  results  can  be  characterized  by. 


E 


SSa 


=  w  E 


SX,m  ^  PiEms)  =  Kx,m  EsX,m-  (76) 

l/^m  I  O 


The  terms  Kx,m  thus  represent  an  integration  efficiency  over  the  input  face, 

Kx,m  =  y"  M  ^  PiEms).  (77) 


l/^m  I  8 


Similar  coefficients  will  exist  for  Y-type  and  Z-type  surfaces.  Optimally,  one 
would  expect  the  Kx,m^  Ky^m^  and  Kz,m  terms  all  to  equal  unity,  but  due  to  the 
choice  of  the  symmetry  set  and  peculiarities  of  the  phase  function  representation 
itself,  suboptimal  conditions  may  exist.  In  our  example,  then,  the  total  flux 
dispersed  at  the  X-type  interface  is  modeled  as  Kx,m  Esx,m- 


Returning  now  to  the  calculation  of  flux  unaccounted  for  due  to  volumetric 
effects,  let  us  generalize  equations  such  as  Eqs.  (71-72)  and  (61),  where  we  let  (for 
example)  1  —  Ax,m  be  the  fractional  flux  on  input  face  type  X  of  stream  m  that 
must  be  accounted  for  via  surface  effects.  Comparison  between  the  development 
above  and  results  such  as  Eq.  (71)  shows  we  must  have  dx,m  Kx,m  =  1  ~  Ax,m] 
that  is. 


dx,m 


1  Ax,m 

E  X ,m 


(78) 
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The  net  result  is  that  at  a  given  output  surface  and  for  a  particular  output 
stream,  the  total  surface  scattering  contribution  to  that  stream  is  determined 
by  summing  over  all  streams  at  that  surface.  The  fractional  contribution  of  each 
stream  will  depend  on  the  Ax,m  value  for  each  input  stream.  Continuing  the 
X-type  surface  example  we  have, 

^Is,X,s  —  ^  ^  '^m  ^  dx,m  Pm{ Ix,m- 

m 


This  surface  incremental  portion  is  then  added  to  the  transmitted  and  volumetric 
scattered  flux  computed  via  the  methodology  described  previously.  Note  that 
ru  is  considered  a  function  of  the  stream  (as  are  d  and  P)  because  the  behavior 
of  a  given  stream  at  a  given  surface  depends  on  the  direction  in  which  the 
input  stream  is  flowing  across  the  surface  (which  volume  element  the  stream  is 
entering),  regardless  of  the  directionality  of  the  output  stream. 

2.6  Direct  Radiation  Considerations 

The  previous  section  assumed  the  radiance  of  diffuse  streams  was  constant 
over  any  input  face  of  interest.  This  assumption  allowed  the  ‘simple’  forms 
for  transmission,  emittance,  and  surface  scattering  fraction  to  be  derived 
analytically  in  the  cases  studied. 

However,  this  assumption  may  not  be  reasonable  for  direct  (solar/lunar) 
radiation,  due  to  the  relative  dominance  of  direct  sources  of  energy  at  visible 
wavelengths.  Typical  boundary  conditions  for  a  4-km  high  atmospheric  volume 
yield  direct  spectral  irradiances  on  the  order  of  1200  to  1400  W/nP-jim  in  the 
visible,  while  the  diffuse  spectral  radiances  are  on  the  order  of  30  to  100  W/m^- 
/im-sr.  The  differences  in  importance  between  these  two  sources  indicates  a 
distinction  should  be  made  in  processing  these  for  radiative  transfer  calculations. 
For  example,  Zardecki  (1995)  proposed  a  separate  LOS  computation  to  trace 
the  optical  depth  of  direct  illumination  reaching  the  center  of  each  cell  in  his 
radiative  transfer  algorithm. 

2.6.1  Cell  Face  and  Volume- Averaged  Illuminations 

In  the  surface-based  model  discussed  here,  it  would  not  be  appropriate  to  trace 
lines  of  sight  to  a  cell  center  but  rather  to  evaluate  the  average  illumination 
on  each  of  the  three  direct  radiation  input  faces  to  each  cell.  This  can  be 
accomplished  by  sampling  over  each  of  these  surfaces  and  forming  an  average 
transmittance  between  the  modeled  volume  upper  boundary  (assuming  the  Sun 
or  Moon  are  the  only  direct  sources  considered)  and  the  specific  cell  wall.  In 
addition,  one  can  determine  the  mean  illumination  provided  to  the  cell  by 
knowing  the  path  length  that  each  sampled  pencil  of  light  traces  out  as  it  passes 
through  the  cell.  Integrating  over  each  input  surface,  one  can  determine  the 
efficiency  of  direct  radiation  entering  each  face  of  each  volume  element.  This 
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indicates  the  fraction  of  direct  energy  scattered  at  each  surface.  From  these 
calculations,  one  may  determine  the  net  direct  surface  scattering  source  term  at 
each  face  within  the  volume. 

Let  the  azimuth  and  zenith  angle  to  the  source  of  direct  radiation  be  specified 
by  (po  and  respectively.  Then  the  direction  of  propagation  of  the  direct 
illumination  is 

Oo  =  {- sin(^o)  sin(6'o),  -  cos(^o)  sin(6'o),  -cos(6'o)}.  (80) 

Here  note  that  the  azimuth  angle  is  given  in  geographical  units  measured  from 
the  North  (the  Y-axis)  and  increasing  toward  the  East  (X-axis). 

For  each  cell,  there  will  be  at  most  three  surfaces  exposed  to  direct  radiation 
(we  ignore  the  finite  solid  angle  subtended  by  the  source  and  treat  it  as  plane 
parallel  energy).  A  series  of  ray  traces  from  different  sample  points  on  the 
surface  upward  to  the  volume  boundary  are  made  for  each  of  the  three  faces  to 
determine  the  average  transmittance  of  direct  radiation  through  the  volume  to 
these  surfaces. 

For  each  sample  point  on  a  cell,  z,  there  will  be  a  net  transmittance  starting 
at  the  overall  volume’s  upper  boundary  and  continuing  through  the  volume  to 
the  cell  wall  of  interest.  Call  this  transmittance  factor  Secondly,  there 

will  also  be  a  distance  that  unscattered  energy  arriving  at  this  sampling  point 
would  travel  to  reach  an  exit  side  of  the  cell  volume  in  question  (Di).  If  the 
sample  volume  has  an  extinction  coefficient  cr,  then  the  direct  radiation  exiting 
the  volume  through  this  pencil  of  light  will  experience  a  total  transmittance  of. 

Tout, I  =  T,n,texp{-aD,).  (81) 

The  mean  illumination  provided  to  the  volume  by  this  pencil  of  light  will  depend 
on  the  number  of  samples  N  taken  over  the  input  face,  the  face  area  A^,  and 
the  cosine  of  the  pencil  with  the  surface  normal  (call  this  |/io|),  Tin,i^  and  Tout,i- 

The  volume  illuminated  by  the  pencil  {Vi)  will  be 

Id  =  A|/io|AViV.  (82) 

The  average  value  of  the  illumination  due  to  this  pencil  of  illumination  will  then 
be 


—  _  ^0  [Tin,i  —  Tout,i]  -Di|/io|A^/A  _  Aq  [Tin,i  —  Tout,i]  |/^o| 

A^  ^A  W’  ^  ^ 

where  the  illumination  over  each  sample  is  averaged  over  the  full  volume  (A^), 
and  where  Aq  is  the  boundary  irradiance  (W/nV-jim)  which  has  an  implicit 
multiplication  by  |i^o|  (the  cosine  of  Oq  with  the  vertical).  The  total  illumination 
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provided  to  the  cell  via  this  input  surface  is  then  found  by  summing  over  all 
sample  illumination  components: 


N 

Ao  =  Y,^o,^.  (84) 

i=l 

In  the  computer  implementation  of  this  approximation,  computations  over  each 
input  face  are  made  at  successively  higher  N  values  until  the  solution  settles  to 
a  stable  result.  The  results  of  these  calculations  are  the  mean  illumination  to 
each  cell  due  to  each  input  surface  (Aq),  the  mean  transmittance  to  each  cell 
surface, 

1  ^ 

Tin  =  (85) 

i=l 

and  the  mean  transmittance  to  the  far  side  of  the  cell  due  to  energy  passing 
through  each  input  surface, 

1  ^ 

T  out  jy  ^  ^ 

'  i=l 

Given  that  the  direct  radiation  is  modeled  as  an  irradiance  (energy  per  square 
area  perpendicular  to  the  vertical  vector,  which  we  then  transform  into  a  result 
appropriate  to  a  plane  perpendicular  to  the  direction  of  propagation  by  dividing 
by  1^0 )  with  a  directionality  in  the  form  of  a  delta  function  in  the  direction  of 
propagation,  then  the  scattering  into  some  output  stream  is  proportional  to 
P(Oj,  Oo)-  Also,  similar  to  the  results  obtained  for  the  diffuse  method,  one 
may  integrate  over  all  output  streams  over  every  face  of  the  volume  element  to 
obtain  the  net  energy  scattered  due  to  input  direct  radiation  entering  along  a 
given  input  face.  The  mathematics  of  this  process  are  similar  to  those  of  the 
diffuse  case  and  will  not  be  repeated  here,  but  the  results  are  similar.  Each 
face  can  be  assigned  a  particular  efficiency  factor  (Aq)  indicating  the  amount  of 
energy  that  must  be  scattered  via  a  surface  process  proportional  to  (1  —  Aq).  We 
thus  have  an  accounting  method  for  volumetric  and  surface-based  contributions 
of  direct  radiation  to  the  diffuse  streaming  radiances. 

2.6.2  Earth  Curvature  Considerations 

One  consideration  of  particular  significance  to  direct  radiation  calculations  is 
how  to  account  for  the  Earth’s  curvature.  Eor  highly  oblique  angles  of  incidence 
of  direct  light  (near  Sunset /Moonset  or  Sunrise/Moonrise),  the  effects  of  the 
curved  Earth  become  significant.  How  does  one  trace  a  path  through  the 
volume  of  near-infinite  length  when  the  Sun  is  at  nearly  90°  zenith  angle? 
In  particular,  it  is  advantageous  when  running  the  radiative  transfer  model  to 
impose  horizontal  periodicity  on  the  modeled  volume.  This  means  potentially 


(86) 


47 


tracing  through  many  copies  of  the  same  volume  in  tracing  a  line  of  sight  toward 
the  Sun. 

To  avoid  the  possibility  of  infinite  loops,  two  steps  are  taken  in  the  computer 
model.  First,  the  code  is  only  callable  from  its  front  end  processing  program  if 
the  Sun  or  Moon  is  above  the  horizon.  Second,  a  correction  is  made  to  account 
for  the  curvature  of  the  Earth  when  ray  tracing.  This  correction  is  based  on 
a  parabolic  approximation  of  the  earth’s  shape  where  it  is  assumed  that  the 
modeled  region  dimensions  are  small  with  respect  to  the  radius  of  the  Earth. 

Eor  this  development,  we  let  the  nominal  radius  of  a  spherical  Earth  be  fixed  at 
Re  =  6371  km.  To  determine  the  amount  of  correction  to  be  made  at  each  step 
in  the  ray  tracing  process,  one  can  imagine  an  observer  at  height  2;  above  sea 
level  observing  the  Sun  in  an  initial  direction  — Oq  =  {“ho,  “ho,  “Co}-  This 
scenario  is  shown  in  the  diagram  in  figure  11. 


Eigure  11.  Modification  of  the  solar  zenith  angle  with  respect  to  the  gravity 
vector  upon  stepping  along  a  path  out  of  the  atmosphere  in  the  direction  of  the 
direct  source. 


Let  us  now  place  an  observer  a  distance  2;  above  the  Earth’s  surface  and  let  us 
define  an  Earth  centered  coordinate  system;  call  it  Xq.  We  can  then  always 
draw  the  x  and  y  axes  such  that  the  observer  position  is  at  the  x-y  origin:  (0, 
0,  2;  +  Re)-  We  now  step  a  distance  D  along  a  line  of  sight  in  the  direction  of 
the  Sun  to  point  (— /ioT>,  —r/oD,  —^qD  +  2;  +  Re)-  Now,  define  the  horizontal 
distance,  H ,  moved  across  the  Earth’s  surface  during  this  step  along  the  line  of 
sight: 

H  =  +  r]^. 
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The  new  location  will  be  a  fraction  of  an  Earth  radius  away  from  the  original 
point  and  thus  will  have  a  different  direction  for  its  gravity  vector. 


We  have  already  seen  the  relationship  between  Oq  and  the  zenith  angle  ^o-  At 
this  new  point  we  can  again  place  an  observer  a  height  z'  above  the  surface  in 
a  new  Earth  centered  coordinate  system  X' .  To  first  order  this  new  height  will 
be 

/  =  z  (87) 


But,  due  to  the  shift  in  position  across  the  Earth’s  surface,  and  the  consequent 
change  in  the  direction  of  the  gravity  vector,  the  zenith  angle  to  the  Sun  will  be 
modified  to 


e',  =  e^-  H/iz  +  Re). 


Therefore,  the  apparent  elevation  angle  of  the  Sun  must  be  modified  continuously 
along  the  line  of  sight  to  simulate  this  change  in  direction  with  respect  to  the 
gravity  vector  at  each  point.  Here  we  assume  that  over  the  span  of  the  particular 
scattering  volume  of  interest  the  change  in  the  zenith  angle  is  so  small  that  it 
can  be  neglected  in  terms  of  propagation  code  accuracy,  and  that  we  may  begin 
each  trace  using  the  same  angle. 

The  scattering  volumes  to  be  studied  are  actually  quite  small  with  respect  to 
this  curvature,  so  the  effects  of  this  correction  are  only  truly  important  when 
computing  for  a  scattering  model  employing  horizontally  periodic  boundary 
conditions.  Eor  example,  consider  a  scattering  volume  that  is  4  km  on  edge 
horizontally  and  4  km  thick  vertically.  Assume  an  observer  at  zero  height  views 
the  Sun  at  the  horizon  (arguably  a  worst  case).  Then,  for  each  pass  through  the 
scattering  volume  we  would  have. 


d6*o  ~ 


4 

6371 


0.6278  mrad. 


This  result,  in  milliradians,  may  appear  very  small,  yet  the  cumulative  effect  is 
such  that  the  line  of  sight  is  5  km  above  the  surface  within  250  km,  and  20  km 
above  the  surface  within  500  km.  This  indicates  that  we  have  bounded  our 
problem,  but  improvements  can  be  made. 


2.7  Output  Face  Calculations 

The  expressions  obtained  thus  far  exhibit  energy  conservation.  They  also  exhibit 
what  is  believed  to  be  the  proper  behavior  of  the  scattered  energy  in  volumetric 
or  surface-based  terms  in  that  the  modeled  energy  is  distributed  according  to 
radiance  by  direction  such  that  it  is  proportional  to  the  phase  function.  But,  it 
remains  to  cast  these  results  as  they  relate  to  the  output  face  diffuse  radiances. 

To  proceed,  we  formulate  an  equation  set  that  accounts  for  all  the  energy  that 
should  contribute  to  each  output  stream.  The  components  are:  (1)  diffuse 
transmitted  energy  from  each  of  the  three  appropriate  input  faces,  (2)  volumetric 
diffuse  scattering  based  on  average  volume  radiance  for  each  stream,  (3) 
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volumetric  black  body  emittance,  (4)  surface-based  scattering  of  energy  from 
all  diffuse  streams  arriving  at  the  given  face,  (5)  volumetric  scattering  of  direct 
irradiance  energy  based  on  volume  averaged  irradiance,  and  (6)  surface-based 
scattering  into  the  given  stream.  This  step  is  termed  the  iteration  on  the 
scattering  source  (Zardecki  (1995);  Lewis  and  Miller  (1984)),  where  at  the  nth 
stage  of  the  process  the  diffuse  energy  incident  at  a  cell’s  boundaries  may  be 
viewed  as  the  total  diffuse  radiance  resulting  from  (n  —  1)  scatterings. 

For  example,  if  we  consider  an  output  X-type  face  for  stream  s  propagating  into 
the  first  octant,  we  have. 


Iout,X,s  =  lTout,X,s  +  (1  “  Tx,s)  La  +  6Is,X,s^ 


(88) 


where  lTout,x,s  is  the  transmitted  radiance, 

^Taut,s  =  Wx,X,sT-X,  +  X  Itn,X,s 
+  Wz,X,s  T-Z,  +  X  Iin,Z,s 

+  Wy,X,sT-Y,  +  X  Itn,Y,s-  (89) 

Here  the  IF’s  are  areal  weighting  factors  introduced  in  table  3.  The  term 
^Is,x,s  represents  the  surface  scattering  contributions  of  both  direct  and  diffuse 
scattering  for  an  X-type  output  surface  for  stream  s. 

And  lastly,  the  limiting  path  radiance  Lg  is 

Lg  =  ZU  —  ^  ^  Wm  P{jljns)Im  B  W  P(flQa)  (AqX  +  AqY  -|-  Aqz),  (90) 


where  =  Aq-Qs  for  the  direct  radiation  scattering,  and  where 

the  A  terms  are  the  mean  volume  illuminations  due  to  direct  energy  entering 
through  the  x,  y,  and  2;  input  face  for  direct  radiation.  Similar  expressions  apply 
for  streams  travelling  in  other  octants  and  out  other  types  of  faces. 

In  conclusion,  the  analysis  leading  to  Eq.  (88)  began  at  the  end  of  section  2.5.4, 
where  we  were  looking  simply  for  an  expression  for  the  limiting  path  radiance, 
Lg.  But  ,  we  found  that  the  volumetric  method  alone  could  not  account  for 
all  the  energy.  Further  consideration  of  the  problem  led  to  a  surface-based 
correction  and  a  means  of  treating  direct  radiance  effects.  We  thus  now 
have  a  full  theoretical  description  of  the  scattering,  emission,  and  absorption 
process.  We  then  conclude  this  chapter  with  a  brief  description  of  the  computer 
implementation  of  this  theory. 
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2.8  The  Atmospheric  Illumination  Module 

The  Atmospheric  Illumination  Module  (AIM)  is  a  computer  implementation 
that  runs  the  radiative  transfer  model  described  in  the  previous  sections  to 
characterize  a  portion  of  the  lower  atmosphere.  AIM  is  thus  an  interface  that 
runs  the  RT  code  at  various  wavelengths,  sets  up  the  input,  and  postprocesses 
the  code  output.  The  actual  RT  code  which  was  described  previously  has 
been  named  BLITS  (Boundary  Layer  Illumination  and  Transmission  Simulation). 
The  AIM  routine  is  designed  to  run  BLITS  in  several  different  modes.  These 
modes  primarily  are  designed  to  speed  its  execution,  but  also  to  permit  greater 
flexibility  in  the  setting  up  of  run  conditions  and  scenarios  that  can  be  processed. 

In  the  immediate  predecessor  to  the  AIM/BLITS  modeling  system,  Zardecki’s 
BLIRB  model  (Zardecki  1995)  combined  some  features  of  the  AIM  preprocessing 
functions  in  with  the  radiative  transfer  modeling  stages.  That  is,  BLIRB  mixed 
aerosol  and  Rayleigh  scattering  properties  within  the  model  itself.  Thus,  when 
aerosol  properties  were  specified,  they  had  to  be  separately  mixed  for  each 
vertically  different  layer  since  the  molecular  properties  varied  with  height.  Each 
cell  was  then  characterized  by  a  given  scatterer  type.  This  allowed  indexing  of 
scattering  properties  by  cell  position. 

This  method  had  to  be  improved  because  BLITS  requires  more  information  to 
run  than  BLIRB.  Notice  that  the  low-density  transfer  Eqs.  (26)  and  (30)  require 
only  values  of  s  and  a.  And,  since  s  is  a  function  of  the  stream  direction, 
it  is  simply  computed  and  stored,  and  a  is  also  a  single  value  per  cell.  But, 
the  BLITS  approach  requires  storage  of  information  on  the  surface  scattering 
efficiencies  (d),  the  face-to-face  and  face- averaged  transmission  factors  (various 
T  terms),  and  the  volumetric  scattering  source  factors  (,  factors). 

By  moving  the  cell  characterization  processing  out  of  the  RT  model  and  into 
a  preprocessing  stage,  it  was  possible  to  improve  the  characterization  of  the 
model  volume  with  fewer  scattering  types.  To  achieve  this  preprocessing  step, 
we  began  from  a  Legendre  expansion  point  of  view.  To  completely  describe  a 
scattering/absorbing/emitting  cell,  all  the  properties  of  the  cell  are  completely 
described  by  the  quantities:  a  and  waXi^  £  =  0  ...  L.  (Xi  is  defined  in  Eq.  (19), 
as  modified  by  the  scaling  transformation.  In  general  L  is  set  to  iV,  the  order 
of  the  stream  expansion.)  Thus,  each  cell  is  completely  described  (for  a  given 
wavelength)  by  a  vector  of  dimension  L  2. 

Using  this  3D  set  of  vectors,  the  main  AIM  processor  runs  a  pattern  recognition 
clustering  model  that  groups  these  vectors  into  ‘classes’  of  vectors,  each 
describable  by  a  mean  vector.  Thus,  AIM  can  transmit  the  mean  vector 
information  for  each  class  to  the  BLITS  code  along  with  an  index  value  which 
assigns  each  cell  to  one  of  the  classes.  The  scattering  properties  for  a  given  cell 
are  then  indexed  to  the  mean  vector  associated  with  that  class.  To  perform 
the  processing,  initially,  all  the  vectors  are  grouped  into  a  single  class.  The 
clustering  algorithm  then  iterates  such  that  in  each  pass  the  class  with  the 
greatest  variance  in  its  components  is  divided  along  the  major  axis  of  its  greatest 
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variance  components,  and  two  daughter  classes  are  produced.  Each  sample  is 
then  compared  against  the  centroids  of  each  class,  and  a  new  class  is  assigned 
based  on  the  distance  of  that  sample  to  the  closest  class  centroid.  Once  all 
samples  have  been  reassigned,  new  centroids  and  variances  for  each  class  can 
be  recomputed  based  on  the  current  elements  of  each  class.  In  each  iteration,  a 
new  class  of  aerosol  is  created.  That  is,  a  new  centroid  vector  describing  aerosol 
properties  in  its  vicinity  is  created.  This  iteration  process  continues  until  a  user 
selected  maximum  number  of  classes  (centroid  vectors)  is  produced,  or  until 
there  is  no  more  variance  in  any  of  the  classes.  The  class  identity  of  each  cell 
at  the  end  of  this  procedure  is  then  passed  to  BLITS  in  order  to  populate  the 
scattering  volume.  The  mean  vector  component  values  for  the  class  centroids 
are  also  passed  to  characterize  each  scattering  aerosol  indexed  in  the  volume 
population  table. 

The  pattern  recognition  preprocessing  stage  of  AIM  thus  optimizes  the 
representation  of  the  scattering  volume  with  respect  to  a  series  of  aerosol  classes. 
These  classes  allow  for  a  reduction  in  the  total  amount  of  data  that  must  be 
stored  by  limiting  the  transmittance,  illumination  coefficients,  etc.  to  only  those 
listed  according  to  aerosol  classes.  However,  additional  data  storage  compression 
is  also  possible  by  exploiting  symmetries  in  each  class  of  stored  aerosol  data. 
These  symmetries  arise  because  each  octant  of  the  angular  scattering  space  will 
have  identical  transmittance  properties  between  input  and  output  faces.  The 
difference  is  in  the  orientation  of  the  input  and  output  faces.  Similarly,  within  a 
single  octant,  due  to  the  discrete  ordinate  technique,  the  diffuse  streams  oriented 
in  primarily  x,  y,  and  2;  directions  will  exhibit  similar  characteristics  for  diffuse 
scattering  and  transmittance  properties. 

To  exploit  these  properties,  all  internal  arrays  are  computed  in  terms  of  the 
so-called  standard  conditions.  This  means  that  all  results  are  computed  for  the 
first  octant  (stream  flowing  in  the  positive  x,  y,  and  2;  directions)  and  oriented 
such  that  the  2;-axis  component,  1^,  is  greater  than  or  equal  to  either  the  x—  or 
y-axis  components.  Results  obtained  are  then  accessed  by  indirect  addressing. 

Though  this  system  only  exhibits  a  maximum  compression  factor  of  24  for  the 
representation  of  cell  scattering  information,  it  does  represent  some  savings. 
However,  by  far  the  largest  contribution  to  the  data  requirements  for  the  model 
are  for  holding  the  positional/directional  data.  Typical  runs  on  a  64  Mb  work 
station  entail  a  32x32x16  x-y-z  cell  structure  with  24-stream  angular  resolution 
and  a  factor  3  multiplication  to  account  for  stream  values  at  each  X,  Y,  and  Z 
input  walls. 
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3.  Monte  Carlo  Model  Comparisons 

The  purpose  of  any  RT  Monte  Carlo  model  is  to  generate  high-accuracy  statistics 
concerning  scattering,  since  the  ability  to  control  the  scattering  effects  is  greater 
than  that  possible  by  solving  the  radiative  transfer  function  directly.  That  is,  the 
ability  to  describe  the  phase  function  is  not  limited  according  to  the  number  of 
streams  that  must  be  represented  in  the  results  if  we  are  only  concerned  with  flux 
outputs.  Therefore,  an  In-house  Monte  Carlo  (IHMC)  code  has  been  developed 
to  characterize  scattering  scenarios  under  various  propagation  conditions. 

The  primary  difficulty  with  a  Monte  Carlo  approach  is  the  time  required  to 
generate  these  statistics  to  a  desired  degree  of  accuracy.  Since  Monte  Carlo 
methods  depend  on  generating  statistics  which  describe  the  scattering  process, 
photon  noise  is  the  controlling  factor  in  determining  the  number  of  iterations 
necessary  to  produce  a  statistical  solution  to  a  specified  degree  of  accuracy.  But, 
this  implies  that  we  would  like  to  know  how  accurate  a  statistic  is  at  any  point 
in  the  solution  process  so  that  we  know  when  to  stop  iterating.  We  therefore 
present  a  method  of  computing  the  solution  accuracy  in  section  3.1.  Section  3.2 
then  provides  a  short  description  of  the  model,  and  section  3.3  discusses  model 
validation  tests. 

3.1  Assessing  Photon  Noise  Statistics 

Due  to  the  weak  law  of  large  numbers  (c.f.  Stark  and  Woods  1986),  the  solution 
with  a  given  number  of  samples  can  be  tested  to  determine  whether  it  exceeds  a 
fixed  degree  of  accuracy.  With  this  knowledge,  a  Monte  Carlo  code  can  be  easily 
written  that  produces  resulting  statistics  to  arbitrary  accuracy.  Of  course,  higher 
degrees  of  accuracy  require  longer  compute  times,  proportional  to  the  number 
of  samples.  Highly  scattering  {zu  sa  1)  dense  media  require  increased  compute 
times  because  the  number  of  collisions  per  photon  prior  to  exiting  the  volume 
is  a  function  of  the  optical  depth.  Depending  on  the  medium,  the  number  of 
collisions  may  be  on  the  order  of  hundreds  per  photon  for  cloud-type  aerosols. 

The  weak  law  of  large  numbers  states  that  for  a  set  of  independent  and  identically 
distributed  (HD)  random  variables  (RVs)  with  mean  /i  and  variance  cr^,  if 
we  do  not  know  the  mean  of  this  distribution  (/i)  a  priori,  we  can  know  that  an 
estimate  of  the  mean,  /ijv  =  X^/N,  will  approach  the  true  mean  as  the 

number  of  samples  (iV)  approaches  infinity.  In  our  Monte  Carlo  propagation 
model,  the  variables  represent  sample  fluxes  exiting  through  the  six  walls  of 
the  sample  volume.  These  X^’s  are  created  by  producing  packets  of  photons, 
1000  in  each  group.  For  each  packet  of  1000  photons,  we  determine  how  many 
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exit  each  wall  and  divide  by  the  total  to  produce  six  statistics.  Each  of  these 
six  is  a  different  random  variable,  but  multiple  packets  should  have  IID  RVs 
for  each  exit  wall.  Since  the  mean  fluxes  from  each  exit  wall  is  what  we  are 
attempting  to  find,  we  can  produce  a  sample  mean  for  each  wall  by  averaging 
the  results  of  multiple  packets  as  above. 

The  next  step  in  the  analysis  involves  the  use  of  the  Chebyshev  inequality.  Here, 
we  treat  /ijv  itself  as  a  random  variable.  According  to  the  weak  law  of  large 
numbers,  the  mean  of  /ijv  should  be  just  ji,  and  the  variance  is  /N .  The 

Chebyshev  inequality  is  then  invoked,  which  furnishes  a  bound  on  the  probability 
of  a  random  variable  deviating  from  its  mean  by  an  amount  8\ 

(90) 

The  term  / 6‘^  thus  furnishes  an  upper  bound  on  the  extent  of  the  error  in 
estimating  the  mean.  If  we  then  fix  the  accuracy  (A^)  that  we  desire  in  the 
answer,  and  evaluate  8  relative  to  the  value  of  the  sample  mean  itself  {8  =  e  /ijv, 
where  e  is  a  fractional  accuracy  required,  assuming  /ijv  >  0),  then  we  know  that 
we  can  ensure  a  given  level  of  accuracy  in  our  solution  by  requiring: 

2 

-  ,,|  >  i)  <  <  A\  (91) 

iV  e  /ijY 

We  know,  then,  that  we  can  stop  acquiring  additional  samples  whenever 

ViVe/iwA  - 

In  the  algorithm  developed,  e  was  set  to  0.02  and  A  was  set  to  0.5.  Due  to  the 
nature  of  the  criteria  above,  it  is  the  product  of  these  two  that  is  significant. 
Thus,  instead  of  considering  the  results  as  a  2-percent  accuracy  (lOOe)  to  a  75- 
percent  confidence  level  (100(1  —  A^)),  the  results  could  equally  be  viewed  as  a 
5-percent  accuracy  to  a  96-percent  confidence  level  (e  =  0.05,  A  =  0.2),  etc.  In 
other  words,  the  solutions  we  seek  have  a  probability  of  less  than  25  percent  of 
deviating  more  than  2  percent  from  the  actual  solution,  and  a  probability  of  less 
than  4  percent  of  deviating  more  than  5  percent  from  the  actual  solution.  In 
the  implementation  these  restrictions  become  even  more  stringent  because  we 
only  consider  the  worst  case  flux.  That  is,  the  smallest  flux  will  usually  have 
the  largest  fractional  error  and  will  require  more  model  iterations  to  satisfy 
the  convergence  criteria  than  fluxes  from  the  remaining  five  faces.  Thus,  the 
accuracy  of  the  remaining  five  faces  will  always  exceed  the  threshold  while  the 
accuracy  of  the  sixth  surface  will  at  least  meet  the  threshold  convergence. 

Note  Anally  that  is  not  known  a  priori  and  is  also  an  estimate,  but  with  the 
number  of  samples  required  under  the  Monte  Carlo  method,  we  generally  have 
a  very  good  estimate  for  this  value  as  well  whenever  we  have  met  the  above 
criterion  for  /ijv  accuracy. 
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3.2  Model  Overview 

The  Monte  Carlo  model  developed  to  implement  the  above  strategy  was  written 
in  the  CH — h  language.  It  was  designed  to  generate  statistical  fractional  energy 
flux  information  for  photons  passing  through  a  scattering  volume  and  out  one 
of  the  six  sides.  CH — h  is  an  object-oriented  programming  language  (c.f.,  Lafore 
1991),  and  as  such,  several  different  object  classes  were  designed  to  simulate 
the  various  phases  of  the  scattering  process.  Several  different  phase  function 
types  were  programmed  into  a  phase  function  object  class.  Options  included  a 
Henyey-Greenstein  phase  function  and  the  Deirmendjian  type-C.l  cloud  phase 
function  for  a  wavelength  of  0.45  /xm.  The  propagation  of  photons  was  handled 
through  the  development  of  a  ray  object  class.  In  this  class,  a  phase  function 
could  be  accessed  to  determine  the  angle  of  scattering  (deflection)  of  the  photon 
from  its  current  path  at  the  point  of  collision.  Other  classes  were  developed  to 
handle  the  bookkeeping  of  photons  passing  through  various  wall  sections  within 
the  scattering  volume,  angular  information,  positional  information,  etc. 

3.3  Model  Validation/Comparisons 

To  verify  the  validity  of  the  IHMC  model,  test  runs  have  been  compared  with  the 
results  of  McKee  and  Cox  (1974).  In  figures  12  through  14,  the  outputs  of  IHMC 
are  compared  with  output  of  the  McKee  and  Cox  (M&C)  Monte  Carlo  model, 
as  derived  from  figures  3  through  5  of  their  paper.  Comparisons  are  made  for 
incident  zenith  angles  of  direct  radiation  of  0°,  30°,  and  60°,  respectively.  The 
original  data  of  M&C  were  modified  to  include  directly  transmitted  radiation, 
based  on  an  estimate  of  the  optical  depth  of  the  cell  for  each  their  data  points. 
These  corrections  were  made  using  the  transmittance  equations  given  in  chapter 
2  and  the  zenith  angle  of  the  direct  radiation  relative  to  the  normal  to  the  top 
of  the  scattering  volume.  From  their  data,  it  appears  that  volume  single-axis 
optical  depths  of  5,  10,  15,  25,  50,  and  73.5  were  used  by  M&C. 

The  M&C  scenarios  consisted  of  a  cubic  volume  filled  with  a  uniform 
conservative  scattering  material  characterized  by  a  phase  function  obtained  from 
Deirmendjian  (1969)  for  a  type  C.l  cloud  at  a  wavelength  of  0.45  /xm.  This  is  a 
phase  function  simulating  a  cumulus-type  water  cloud.  Conservative  scattering 
was  assumed  in  all  calculations. 

Comparison  of  the  results  in  the  three  figures  indicates  IHMC  closely  follows  the 
M&C  data.  Slight  differences  are  seen  for  the  energy  exiting  the  volume  sides 
for  the  60°  zenith  angle  case,  but  here  the  difference  may  be  in  the  corrections 
applied  to  account  for  the  amount  of  direct  energy  passing  through  the  volume. 
The  exact  values  of  the  optical  depths  used  by  M&C  were  not  reported,  except 
for  the  maximum  optical  depth  case  (73.5).  Another  explanation  is  that  the 
M&C  data  may  be  slightly  in  error:  Davies  (1978)  reported  small  discrepancies 
when  he  compared  his  own  results  with  the  M&C  calculations  and  attributed 
the  differences  to  weaknesses  in  the  M&C  algorithm.  Note,  for  example,  that  in 
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Figure  12.  Scattered  and  transmitted  energy  exiting  from  a  cubical  volume. 
Symbols  represent  M&C  data  modified  to  include  transmitted  direct  energy. 
Curves  represent  IHMC  results.  Incident  energy  is  a  plane  parallel  beam  striking 
the  top  of  the  volume  with  a  normal  incidence. 


Figure  13.  Monte  Carlo  comparisons  similar  to  previous  figure  except  incident 
energy  strikes  the  top  with  a  zenith  angle  of  30°  relative  to  the  normal.  Incident 
radiation  strikes  top  and  one  side  of  the  volume. 
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figure  14  the  geometric  result  for  energy  exiting  through  the  sides  at  zero  optical 
depth  should  be  63.4  percent.  The  trend  in  the  corrected  M&C  data  is  for  less 
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Figure  14.  Results  similar  to  previous  figure  except  that  incident  energy  strikes 
the  top  with  an  angle  of  60°  relative  to  the  normal.  Incident  radiation  strikes 
top  and  one  side  of  the  volume. 

than  60  percent  to  be  flowing  out  the  sides.  Similarly  the  geometric  result  for 
energy  exiting  the  volume  base  should  be  36.6  percent.  Both  of  these  targets 
appear  consistent  with  the  trends  at  low  optical  depth  for  the  in-house  model. 

From  these  figures,  there  is  a  determination  that  IHMC  achieves  a  reasonable 
performance  when  compared  to  the  M&C  results.  Further,  since  we  can 
control  the  exact  model  specifications  for  the  cases  of  comparison,  there  is  no 
question  regarding  what  optical  depths  to  apply  in  running  the  BLITS  code. 
The  IHMC  model  will  thus  be  used  henceforth  for  comparisons  with  the  BLITS 
code.  However,  we  will  restrict  the  application  of  IHMC  to  only  those  scenarios 
described  by  M&C,  for  which  we  have  direct  verification  that  the  code  is  valid. 
(Note:  we  compared  results  using  the  BLITS  model  because  M&C  dealt  with 
specific  uniform  density  media.  AIM  generates  input  data  for  general  atmospheric 
conditions  and  would  be  inappropriate  for  comparison.) 

3.4  Analysis 

The  results  of  the  BLITS  code  were  compared  to  those  of  the  IHMC  code  for 
each  scenario  studied  by  M&C.  Results  were  obtained  from  the  Monte  Carlo 
code  using  the  full  Diermendjian  C.l  phase  function.  BLITS  results  were 
produced  using  modified  Legendre  expansion  coefficients  of  the  phase  function 
such  that  the  forward  peak  was  removed  via  one  of  two  different  scaling 
transformations.  The  first  transformation  was  the  (5-M  method  proposed  by 
Wiscombe  (1977).  The  second  was  the  log-least-squares  (LLS)  method  described 
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in  the  appendix.  The  second  method  was  developed  to  suppress  higher  order 
Legendre  components,  thereby  avoiding  negative  values  in  the  phase  function 
approximation.  These  negative  values  can  lead  to  erroneous  predictions  of 
stream  radiances. 

Comparisons  were  made  for  various  levels  of  resolution  of  the  calculation  grid 
for  the  cubic  scattering  volumes  used  by  M&C.  These  calculations  show  the 
overall  robustness  of  the  method.  The  walls  of  the  volume  were  considered  to 
be  bounded  by  a  non-reflecting/non-emitting  void.  The  external  illumination 
consisted  of  a  plane  parallel  source  incident  on  the  top  and  one  side  of  the 
volume  (virtually  a  cubical  gas  in  deep  space).  The  cases  considered  varied 
the  medium’s  single-axis  optical  depth  and  the  zenith  angle  of  arriving  direct 
radiation.  Results  were  computed  for  direct  energy  incident  at  the  top  of  the 
volume  and  striking  the  eastern  side  from  0°,  30°,  and  60°  zenith  angles.  We 
did  not  use  the  M&C  data,  except  for  comparison.  Instead,  we  computed  results 
using  IHMC  and  BLITS  for  single  axis  optical  depths  of  the  medium  of  5,  7.5,  10, 
12.5,  15,  20,  25,  37.5,  50,  and  75. 

In  figure  15  we  compare  the  0°  zenith  angle  incidence  direct  source  case  for  flux 
reflected  off  the  top  of  the  modeled  volume  using  IHMC  and  BLITS  with  varying 
numbers  of  cells.  This  figure  is  based  on  BLITS  runs  using  the  (5-M  method, 
and  the  results  show  that,  even  with  a  nominal  35  optical  depths  per  cell,  the 
2x2x2  coarseness  calculations  produce  reasonable  results.  Figure  16  shows  the 
same  output  for  the  LLS  method  described  in  the  appendix.  Note  that  the  LLS 
method  appears  not  to  predict  results  as  well  as  (5-M  at  the  2x2x2  resolution, 
but  appears  to  approach  the  IHMC  results  asymptotically  as  the  resolution  is 
increased. 

Figures  17  through  20  show  details  of  the  fluxes  emitted  from  the  sides  and 
bottom  of  the  scattering  region  at  different  computation  resolutions  for  the 
(5-M  and  LLS  methods  for  the  same  0°  incident  scenario.  As  can  be  seen, 
while  the  performance  of  the  algorithm  appears  robust  and  seems  to  operate 
efficiently  for  both  the  (5-M  and  LLS  methods,  the  LLS  approach  appears  to 
slightly  outperform  the  (5-M  results  at  the  higher  resolutions,  which  we  assume 
behave  in  an  asymptotic  fashion.  Similar  results  were  obtained  for  the  M&C 
cases  of  30°  and  60°  angles  of  incidence  of  the  direct  radiation  beam.  In  all  cases 
studied,  the  LLS  method  appeared  to  outperform  the  traditional  (5-M  approach. 
The  primary  factor  for  this  increase  in  performance  is  attributed  to  the  nature  of 
the  phase  function  and  the  limitation  imposed  that  only  a  24-stream  model  was 
used.  For  this  case,  the  (5-M  corrected  Legendre  expansion  of  the  phase  function 
produces  negative  phase  function  results  in  8  out  of  24  different  scattering 
directions.  These  negative  values  produce  negative  stream  values  after  the  direct 
radiation  scatters  into  these  directions.  Once  these  negative  radiances  are  in 
place,  they  then  disrupt  the  overall  evaluated  flux  statistics.  Introducing  more 
streams  would  solve  the  problem  of  the  negative  radiances  because  a  higher 
order  Legendre  expansion  could  be  used,  but  this  would  require  larger  system 
memory  and  slower  model  operation. 
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Figure  15.  Comparisons  for  the  (5-M  method  at  varying  cell  resolutions  of  the 
scattering  volume  for  0°  zenith  angle  incident  direct  energy  for  energy  exiting 
volume  top. 


Figure  16.  Comparisons  for  the  LLS  method  at  varying  cell  resolutions  of  the 
scattering  volume  for  0°  zenith  angle  incident  direct  energy  for  energy  exiting 
volume  top. 
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Figure  17.  Comparisons  for  the  (5-M  method  at  varying  cell  resolutions  of  the 
scattering  volume  for  0°  zenith  angle  incident  direct  energy  for  energy  exiting 
volume  sides. 


Figure  18.  Comparisons  for  the  LLS  method  at  varying  cell  resolutions  of  the 
scattering  volume  for  0°  zenith  angle  incident  direct  energy  for  energy  exiting 
volume  sides. 
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Figure  19.  Comparisons  for  the  (5-M  method  at  varying  cell  resolutions  of  the 
scattering  volume  for  0°  zenith  angle  incident  direct  energy  for  energy  exiting 
volume  base. 
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Figure  20.  Comparisons  for  the  LLS  method  at  varying  cell  resolutions  of  the 
scattering  volume  for  0°  zenith  angle  incident  direct  energy  for  energy  exiting 
volume  base. 


61 


However,  disregarding  for  the  moment  any  differences  between  the  LLS  and  S- 
M  performances,  the  primary  characteristic  of  these  results  is  that  the  BLITS 
scattering  algorithm  appears  to  approximate  the  Monte  Carlo  model  outputs 
to  a  high  accuracy  at  all  optical  depths  for  the  4^  and  16^  cases.  The  worst 
case  (the  0°  flux  exiting  the  volume  sides)  fails  at  the  two-cell  approximation 
simply  because  two  cells  per  axis  are  insufficient.  Certainly  the  four-cell  case 
has  much  better  performance  at  70  optical  depths  than  does  the  two-cell  model 
at  17  optical  depths,  even  though  these  two  cases  have  roughly  the  same  optical 
depth  per  cell.  This  indicates  that  the  overall  success  of  the  model  seems  to 
depend  on  the  feedback  produced  by  sufficient  numbers  of  lattice  walls  interior 
to  the  modeled  volume. 

3.5  3D  Segmented  Cloud  Cases 

Analysis  of  the  BLITS  model  performance  in  the  M&C  cases  revealed  that  two 
cells  per  axis  was  insufficient  to  fully  test  the  viability  of  the  model,  since  the 
model  relies  on  interaction  between  scattering  effects  at  different  interior  walls 
within  the  modeled  volume.  But,  with  the  elimination  of  the  single  cell  and 
two-cell-per-axis  cases,  the  maximum  optical  depth  per  cell  that  is  treated  is 
75/4=18.75.  This  number  is  relatively  low  when  compared  to  the  optical  depth 
of  a  typical  scale  cloud  puff  of  1/4  km  on  a  side.  In  order  to  test  the  model  at 
higher  aerosol  concentrations,  a  set  of  modified  scenarios  had  to  be  developed. 
But,  these  new  scenarios  should  have  certain  characteristics.  It  would  be  possible 
simply  to  run  uniform  volume  cases  like  the  M&C  scenarios,  but  at  much  higher 
optical  depths.  And,  while  this  course  has  some  merit,  it  does  not  address 
the  issue  of  multiple  scattering  interaction  between  different  cloud  elements. 
Further,  there  is  no  consistency  check  on  the  results  obtained  to  ensure  that  there 
is  not  some  flaw  in  the  technique  at  higher  optical  depths.  Another  possibility 
is  to  devise  a  means  of  concentrating  the  aerosols  into  fewer  cells,  creating  gaps 
between  the  different  cloud  elements,  and  thus  permit  testing  of  cloud-to-cloud 
scattering  interactions. 

The  result  of  these  considerations  was  a  method  which  conformed  to  the  M&C 
scenarios  in  one  limit,  but  which  modeled  separate  cloud  sections  as  a  function 
of  a  pair  of  related  separation  parameters,  q  and  P. 

Let  the  overall  modeled  volume  be  described  by  a  cube  of  unit  length  (1  km) 
on  each  axis.  Then,  define  the  quantity  g  as  a  fractional  length  and  divide  the 
volume  into  eight  cubes,  one  each  with  an  outer  corner  positioned  in  each  corner 
of  the  overall  volume,  as  shown  in  figure  21.  The  parameter  q  measures  the 
single-axis  fractional  length  of  each  of  the  eight  scattering  subregions  within  the 
overall  volume  (0  <  g  <  1/2).  The  remainder  of  the  overall  volume  is  considered 
vacuum.  The  fractional  distance  along  a  single  axis  from  center-to-center  of  two 
of  these  cubes  is  1  —  q.  We  then  introduce  the  parameter  P  as  P  =  q/(l  —  g), 
which  varies  between  0,  where  all  the  scattering  material  has  been  compressed 
into  the  eight  corner  regions,  and  1,  where  the  scattering  material  is  uniformly 
distributed  over  the  entire  scattering  volume.  Thus  P  =  1  corresponds  to  the 
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M&C  scenarios  and  P  =  0  corresponds  to  a  vacuum  condition.  We  thus  have 
complete  information  about  the  bounding  behavior  of  the  scenarios  proposed  and 
can  investigate  the  intermediate  behavior  to  see  that  it  smoothly  transitions  by 
varying  the  P  parameter. 


1^ 


_ 


P  =  q/(l-q) 


Figure  21.  Cloud  subregion  configuration  within  the  overall  scattering  volume 
utilizing  q  and  P  parameters. 


Figures  22  through  30  illustrate  the  different  fluxes  (top,  base,  and  sides)  for 
the  three  different  incident  flux  zenith  angles.  The  optical  depths  in  each 
scenario  have  been  normalized  to  the  uniform  density  scatterer  case  such  that  a 
constant  amount  of  scatterer  is  present  in  the  overall  volume  as  P  varies.  This 
is  accomplished  by  computing  a  modified  extinction  coefficient,  a'  =  a/(8q^), 
which  accounts  for  the  concentration  of  the  scattering  materials  in  the  eight 
corner  regions.  Note  that  as  P  decreases  from  unity  more  of  the  incident 
radiation  passes  unscattered  through  the  volume.  This  is  significant  in  the 
analysis  of  results. 

Using  these  Monte  Carlo  results  as  the  baseline,  we  then  proceeded  to  run  the 
BLITS  model  at  a  series  of  resolutions  (number  of  cells  per  axis)  for  each  scenario. 
Similar  to  the  previous  figures,  we  ran  the  model  at  1,  2,  4,  8,  16,  and  32 
cell  resolutions,  hereafter  referred  to  as  the  U,  2^,  4^,  8^,  16^,  and  32^  cases, 
respectively,  in  reference  to  the  total  number  of  cells  modeled  in  the  volume.  We 
used  the  same  average  densities  per  optical  axis  as  in  the  previous  cases  run,  and 
set  the  P  parameter  to  1/3,  3/5,  and  7/9,  corresponding  to  q  values  of  1/4,  3/8, 
and  7/16,  respectively.  Note  that  in  all  three  of  the  new  P  value  cases  (P  =  1  in 
the  original  M&C  scenarios),  the  2^  model  resolution  was  incapable  of  resolving 
the  density  variations.  In  terms  of  the  analysis  below,  these  cases  are  termed 
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Fractional  Total  Energy 


Figure  22.  Fractional  flux  escaping  from  volume  top  for  0°  incident  direct 
radiation  as  a  function  of  the  P  parameter  and  mean  single-axis  optical  depth. 


Fractional  Total  Energy 


Figure  23.  Fractional  flux  escaping  from  volume  sides  for  0°  incident  direct 
radiation  as  a  function  of  the  P  parameter  and  mean  single-axis  optical  depth. 
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Figure  24.  Fractional  flux  escaping  from  volume  base  for  0°  incident  direct 
radiation  as  a  function  of  the  P  parameter  and  mean  single-axis  optical  depth. 
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Figure  25.  Fractional  flux  escaping  from  volume  top  for  30°  incident  direct 
radiation  as  a  function  of  the  P  parameter  and  mean  single-axis  optical  depth. 
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Figure  26.  Fractional  flux  escaping  from  volume  sides  for  30°  incident  direct 
radiation  as  a  function  of  the  P  parameter  and  mean  single-axis  optical  depth. 


Figure  27.  Fractional  flux  escaping  from  volume  base  for  30°  incident  direct 
radiation  as  a  function  of  the  P  parameter  and  mean  single-axis  optical  depth. 
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Figure  28.  Fractional  flux  escaping  from  volume  top  for  60°  incident  direct 
radiation  as  a  function  of  the  P  parameter  and  mean  single-axis  optical  depth. 


Figure  29.  Fractional  flux  escaping  from  volume  sides  for  60°  incident  direct 
radiation  as  a  function  of  the  P  parameter  and  mean  single-axis  optical  depth. 
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Figure  30.  Fractional  flux  escaping  from  volume  base  for  60°  incident  direct 
radiation  as  a  function  of  the  P  parameter  and  mean  single-axis  optical  depth. 


under-resolved.  The  4^  model,  on  the  other  hand,  was  able  to  critically-resolve 
the  P  =  1/3  case.  That  is,  in  the  4^  model  each  cloud  segment  occupied  exactly 
1  cell  in  the  corner  of  the  volume.  In  all  other  cases,  the  4^  model  under-resolved 
the  cloud  segments.  For  the  8^  model,  in  the  P  =  1/3  case,  each  cloud  segment 
occupied  a  2^  cell  region  in  each  corner  of  the  volume,  and  in  the  P  =  3/5  case 
each  cloud  segment  occupied  a  3^  region  in  each  corner.  But,  the  8^  model 
was  unable  to  resolve  (under-resolved)  the  cloud  edges  in  the  P  =  7/9  case. 
We  thus  were  able  to  study  several  conditions:  under-resolved  cases  where  the 
physical  model  was  unable  to  resolve  the  density  variations,  critically-resolved 
cases  where  the  physical  boundaries  of  the  cells  could  just  match  the  boundaries 
of  the  cloud  segments  (for  example,  in  the  P  =  1/3  case  using  the  4^  model),  and 
those  cases  that  were  over-resolved  (for  example,  the  32^  model  for  the  P  =  3/5 
case  where  both  the  regions  containing  the  cloud  segments  and  the  gaps  between 
the  segments  were  able  to  be  modeled  as  more  than  one  cell  wide).  The  behavior 
of  BLITS  under  these  various  conditions  allows  the  testing  of  the  robustness  of 
the  code  in  providing  flux  data  under  various  degrees  of  stressing  conditions. 

However,  following  computation  of  these  scenarios,  it  became  apparent  that 
except  for  the  P  and  2^  models  which  had  already  been  shown  to  lack  sufficient 
numbers  of  cell  walls  to  produce  reliable  results,  the  remainder  of  cases  fell 
decisively  into  only  two  groups:  a  first  group  containing  all  the  under-resolved 
cases,  and  a  second  group  containing  both  the  critically-resolved  and  over¬ 
resolved  cases.  The  under-resolved  cases  consisted  of  the  4^  model  results  for  the 
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P  =  3/5  and  P  =  7/9  cases  and  the  8^  model  results  for  the  P  =  7/9  case.  In  all 
the  critically-resolved  and  over-resolved  cases,  the  model  exceeded  expectations 
in  producing  very  close  results  to  the  Monte  Carlo  run  outputs  for  cell  optical 
depths  up  to  37.5  per  cell  (pre-scale  transform  optical  depth).  Even  after  the 
scale  transform,  which  reduced  the  extinction  coefficient  by  about  60  percent, 
we  were  still  treating  up  to  20  optical  depths  per  cell  in  the  worst  case  (the  4^ 
model  in  the  P  =  1/3  case).  For  this  case,  we  had  an  overall  RMS  error  of 
1.57  percent  of  the  total  flux.  The  results  for  all  cases  are  shown  in  table  5. 

Table  5.  RMS  errors  for  all  cases  run;  under-resolved  in  brackets,  critically- 
resolved  in  parentheses,  over-resolved  without  brackets. 


p 

%  RMSE  1 

R 

2^ 

43 

8^ 

16^ 

32^ 

1/3 

[22.8] 

[25.3] 

(1.67) 

0.76 

0.638 

0.534 

3/5 

[11.4] 

[13.4] 

[12.9] 

(0.80) 

0.714 

0.597 

7/9 

[C.31 

[C.31 

[7.06] 

[6.86] 

(0.901) 

0.858 

1 

(7.7) 

2.5 

0.98 

0.93 

0.857 

0.934 

To  analyze  the  significance  and  meaning  of  these  results,  we  divide  the 
consideration  of  these  results  into  low-,  medium-,  and  high-resolution  cases. 
Figures  31  through  34  show  the  low  resolution  cases  of  the  R  and  2^  models. 
Since  they  are  only  resolved  in  the  uniform  density  (P  =  1)  case,  they  illustrate 
the  situation  of  a  cloud  with  spatial  structure,  which,  due  to  the  model  spatial 
resolution,  nevertheless  occupies  only  a  single  cell  in  the  modeled  space.  One 
interesting  feature  of  the  table  is  a  comparison  between  the  2^  and  the  single 
cell  model  results  for  the  P  =  1/3  and  P  =  3/5  cases.  Here,  the  single  cell 
actually  outperforms  the  2^  model.  The  reason  appears  to  be  that  the  2^  results 
do  a  much  better  job  of  characterizing  the  P  =  1  case  to  the  extent  that  they 
are  over-optimized.  That  is,  in  a  sense,  the  R  model  appears  to  be  doing  the 
job  of  representing  an  ensemble  of  different  density  conditions  that  all  map  onto 
a  uniform  density  problem  when  expressed  in  a  single  cell  model.  Of  course, 
the  low  P  cases  have  the  volume  almost  empty  of  material  and  so  represent 
worst  case  conditions  of  applying  the  uniform  cube  model  as  a  representation  of 
physical  reality. 

The  next  several  cases  consider  the  various  results  of  the  4^  and  8^ 
characterizations.  These  cases  are  interesting  because  they  provide  a  better 
insight  into  the  model  behavior  due  to  their  greater  numbers  of  interior  cell 
walls.  They  also  illustrate  critically  resolved  cases  of  P  =  1/3  using  the  4^ 
model  and  P  =  3/5  using  the  8^  model,  both  of  which  show  excellent  response 
under  high-density  cell  conditions.  These  results  are  given  in  Figures  35  through 
40.  One  of  the  most  interesting  features  is  the  similarity  between  results  in  the 
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Figure  31.  Uniform  density  media  for  the  one  cell  model. 


Figure  32.  Under-resolved  cases  for  the  one  cell  model. 

two  cases,  for  even  though  the  8^  model  has  additional  flexibility,  the  results 
for  the  unresolved  P  =  7/9  case  are  very  similar  in  both  examples  even  though 
the  8^  model  has  a  few  of  its  inner  cells  with  somewhat  lower  optical  depths. 
These  results  tend  to  indicate  that  natural  clouds  will  be  so  optically  thick  that 
their  influence  will  extend  outside  their  volume  if  too  coarse  of  a  grain  is  used  in 
modeling  the  geometry  of  the  cloud  cellular  structure.  These  extensions  might, 
under  unlucky  circumstances,  result  in  inappropriate  fluxes  reaching  the  surface. 
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Figure  33.  Uniform  density  media  case  for  the  2^  cell  model. 


Figure  34.  Under-resolved  cases  for  the  2^  cell  model. 

In  most  cases,  then,  the  effects  of  coarse  grids  will  be  an  increase  in  estimated 
cloud  layer  reflectivity  and  a  decrease  in  estimated  net  layer  transmittance. 
These  observations  have  obvious  implications  beyond  the  domain  of  our  present 
problem  because  they  indicate  some  of  the  weaknesses  of  ID  RT  models,  in  that 
these  models  will  tend  to  underestimate  the  net  solar  loading  to  the  ground  and 
thus  may  influence  global  energy  balance  estimates. 

The  final  group  of  cases  are  the  high  resolution  16^  and  32^  cell  models.  Here, 
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Figure  35.  Critically-resolved  cases  for  the  4^  cell  model. 


Figure  36.  Resolved  cases  for  the  4^  cell  model. 

there  is  only  one  critically-resolved  scenario  (P  =  7/9  for  the  16^  model), 
with  the  remainder  over-resolved.  As  anticipated,  these  results  all  show  a  high 
correlation  to  the  Monte  Carlo  results.  Here,  we  also  note  that  it  would  be 
difficult  to  actually  produce  results  that  are  much  closer  to  the  Monte  Carlo 
model  outputs  regardless  of  the  number  of  cells  we  used  unless  we  were  to 
rerun  the  Monte  Carlo  cases.  Since  we  only  used  2  percent  accuracy  in  the 
Monte  Carlo  runs,  it  is  unlikely  we  could  improve  much  on  the  approximately 
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Figure  37.  Under-resolved  cases  for  the  4^  cell  model. 


Figure  38.  Critically-resolved  cases  for  the  8^  cell  model. 

1  percent  accuracy  of  the  RMSE.  Note  further  that  in  all  the  16^  and  32^  results, 
the  correlation  to  the  Monte  Carlo  results  was  0.999  or  greater,  which  validates 
both  the  Monte  Carlo  approach  followed  as  well  as  the  BLITS  model  results. 
The  results  for  these  comparisons  are  found  in  figures  41  through  43. 
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Figure  39.  Resolved  cases  for  the  8^  cell  model. 


Figure  40.  Under-resolved  case  {P  =  7/9)  for  the  8^  cell  model. 
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Figure  41.  Critically-resolved  case  (P  =  7/9)  for  the  16^  cell  model. 


Figure  42.  Resolved  cases  for  the  16^  cell  model. 
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BLITS  Percent  Flux 
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O  Resolved  P  =  1/3  cases,  RMSE  =  0.53% 
A  Resolved  P  =  3/5  cases,  RMSE  =  0.60% 
^  Resolved  P  =  7/9  cases,  RMSE  =  0.86% 
80  I  Q  Resolved  P  —  1  cases,  RMSE  —  0.93% 


Monte  Carlo  Percent  Flux 


Figure  43.  All  cases  for  the  32^  cell  model. 


4.  Using  AIM  Output 


The  primary  purpose  of  the  RT  codes  discussed  in  this  report  is  to  facilitate 
the  radiometrically  correct  rendering  of  the  appearance  of  natural  cloud  scenes. 
This  task  encompasses  several  related  subtasks.  For  example,  one  may  be 
interested  in  the  appearance  of  a  simulated  surface  below  a  broken  cloud  field. 
One  may  wish  to  render  objects  (trees,  buildings,  bridges,  vehicles)  that  are 
either  attached  to  or  moving  on  or  near  the  terrain.  One  may  wish  to  determine 
the  atmosphere’s  appearance  (fog,  haze,  or  precipitation)  under  clouds.  Or,  one 
may  wish  to  characterize  the  appearance  of  the  cloud  field  itself. 

Each  of  these  subtasks  entails  a  sequence  of  different  processing  steps.  However, 
this  chapter  focuses  on  the  rendering  of  components  of  the  atmosphere.  For 
this  set  of  tasks,  one  must  perform  integrations  over  specified  lines  of  sight. 
Each  LOS  will  have  an  initial  radiance  value  plus  path  integrated  properties. 
These  path  properties  entail  transmittance  and  path  radiance  characteristics. 
The  path  radiance  requires  interpolations  over  direction  and  position  of  RT  code 
results  obtained  in  a  set  of  specific  directions  at  regularized  positions.  Positional 
interpolation  can  be  accomplished  any  number  of  ways,  though  the  process  is 
generally  based  on  distance  to  nearest  neighbor  points.  Directional  interpolation 
depends  on  the  Gaussian  quadrature  method  chosen  to  simulate  the  scattering. 
This  chapter  ends  with  a  description  of  a  cloud  rendering  methodology  and 
several  output  file  formats  useful  for  input  to  visualization  routines. 

4.1  Line  of  Sight  Calculations 

For  any  medium  we  have  a  standard  technique  for  representing  the  effects  of 
the  volume  for  a  given  LOS  through  the  volume.  This  technique  is  related  to 
the  equation  of  radiative  transfer  (22),  and  its  solution,  Eq.  (24).  Since  any 
image  to  be  rendered  can  be  considered  composed  of  a  series  of  pixels,  and  each 
pixel  can  be  thought  of  as  representing  a  particular  LOS  through  the  medium, 
determining  the  atmospheric  effects  for  a  particular  image  becomes  a  matter 
of  tracing  a  large  number  of  lines  of  sight  through  the  simulated  atmosphere. 
Each  LOS  is  characterized  by  an  origin,  a  direction,  and  a  distance,  and  effects 
are  determined  by  tracing  through  a  sequence  of  adjacent  cells  in  the  scattering 
volume. 
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4.2  Integration  Procedure 

For  the  moment,  let  us  assume  that  we  know  the  original  value  for  the  scene 
radiance  at  the  far  end  of  each  LOS  associated  with  each  image  pixel.  All 
that  remains  then  is  to  evaluate  the  transmittance  and  integrated  path  radiance 
associated  with  the  intervening  atmosphere.  Call  the  initial  radiance  associated 
with  a  pixel  Iq.  Let  S  represent  the  path  length  of  the  intervening  atmosphere. 
Let  Ts  be  the  transmittance,  and  Ls  be  the  integrated  path  radiance  of  that 
atmosphere.  Then,  if  Is  is  the  radiance  presented  at  the  terminus  of  that 
atmospheric  path  as  the  energy  enters  the  imaging  optics,  we  have. 

Is  =  IqTs  +  Ls-  (93) 

The  entire  effect  of  the  atmospheric  path  can  thus  be  represented  by  Ts  and  Ls- 

While  some  simplistic  models  compute  Ts  and  Ls  from  single  values  of  extinction 
coefficient  and  limiting  path  radiance  (as  in  the  Introduction),  these  are  actually 
integrated  quantities.  Normally,  atmospheric  path  effects  are  evaluated  by 
dividing  the  LOS  into  a  series  of  small  increments  from  the  observed  object 
to  the  observer.  But  these  equations  can  also  be  evaluated  outward  from  the 
observer,  and  we  shall  show  why  this  method  has  certain  advantages  during  the 
rendering  process. 

First,  we  consider  the  traditional  method  running  the  LOS  from  the  object  to 
the  observer.  Let  n  be  an  index  that  runs  from  1  to  iV,  and  let  Sn  be  the 
length  of  a  given  path  increment  along  a  line  from  the  object  to  the  observer. 
We  then  define  as  the  transmittance  over  the  path  increment  Sn  and  In  as 
the  limiting  path  radiance  over  that  increment.  Since  we  have  now  converted  a 
continuous  calculation  into  a  series  of  increments,  let  us  define  Tjv  =  Ts  as  the 
total  transmittance  over  the  entire  LOS,  and, 

Tn  =  Tn-l  in',  Tq  =  1.  (94) 

Similarly,  let  us  define  Sn  as 

Sn  =  +  ^n',  So  =  0.  (95) 

This  implies  Sn  =  S.  Lastly,  we  have  a  recursive  definition  for  Ln, 

Ln  —  Ln  —  1  tn  T  0  (1  tn'),  Lq  —  0,  Ln  —  Ls-  (^5) 

Here,  To,  5*0,  and  Lq  are  initial  values.  However,  it  is  often  more  valuable  to 
determine  the  path  characteristics  from  the  observer’s  perspective  by  integrating 
the  path  effects  beginning  at  the  observation  point  and  proceeding  out  along  the 
LOS.  For  this  consideration,  let  us  define  a  new  index  variable  m  =  1  +  iV  —  n. 
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We  can  then  define  a  new  set  of  transmittance  and  path  radiance  variables, 
m?  trn  —  ^I+W— m?  and  _l_ jy _ •  Then, 

S.'..  =  5:„-i+4;  S'„  =  0.  (97) 

=  ^m-1  Tq  =  1.  (98) 

L'm  =  -^m-1  +  ^’m-1  C  -^0  =  O'  (99) 

This  second  path  description  has  the  advantage  that  objects  at  different 
ranges  along  the  same  LOS  can  be  interpolated  using  results  at  different  m 
indices.  These  results  thus  would  work  well  for  moving  objects  as  they  change 
ranges  and  possibly  position  within  the  field  of  view,  assuming  a  stationary 
observer.  A  further  advantage  of  this  approach  is  that  if  there  is  some  threshold 
transmittance  below  which  we  are  essentially  viewing  only  atmosphere,  then  the 
path  integration  can  be  truncated,  saving  computer  resources.  In  contrast,  the 
disadvantage  of  the  former  approach  is  that  even  if  the  observer  is  moving,  there 
will  only  be  at  most  a  single  LOS  that  obtains  the  advantage  that  results  can 
be  reused  for  the  observer  at  different  ranges  from  the  object.  All  other  lines 
of  sight  must  be  continually  recomputed.  It  therefore  appears  that  the  second 
approach  is  more  versatile. 

To  evaluate  the  quantities  is  relatively  routine.  Over  each  path  segment,  there 
will  be  a  given  extinction  coefficient,  associated  with  the  cell  in  the  modeled 
domain  being  traversed.  The  transmittance  over  that  path  segment  will  then  be 
=  exp(  —  k'^s'^).  The  distances  traversed  in  each  cell  can  be  determined 
through  a  ray  tracing  algorithm.  For  cubical  cell  shapes,  the  computation 
of  s'^  is  relatively  straightforward.  However,  due  to  the  computer-intensive 
computation  of  it  may  be  more  efficient  to  use  less  exact  means  whereby 
the  path  is  sampled  in  evenly  spaced  increments.  This  latter  method  allows  for 
more  rapid  simulation  of  the  path  geometry  but  could  lead  to  difficulties  when 
the  observer  is  in  motion  with  respect  to  the  cloud  features. 

The  computation  of  the  terms  is  somewhat  more  complicated.  We  treat  this 
in  two  phases.  First,  volumetric  limiting  path  radiances  must  be  derived  from 
the  surface  based  radiance  results  of  the  RT  code.  This  first  stage  produces 
results  at  the  specific  discrete  ordinate  directions  for  which  the  code  was  run. 
These  results  must  then  be  interpolated  by  direction  to  produce  a  limiting  path 
radiance  for  the  specific  direction  leading  into  the  observer’s  sensor.  These  topics 
are  treated  in  the  next  two  sections. 
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4.3  Evaluating  Limiting  Path  Radiance 

The  extraction  of  limiting  path  radiance  from  the  RT  code  results  requires 
some  reorientation  in  viewpoint,  since  the  RT  code  output  is  surface  averaged 
radiances,  while  the  integrated  path  radiance  calculation  requires  volume-based 
information. 

In  order  to  evaluate  this  information,  we  consider  the  energy  of  a  single  stream 
(index  i)  as  it  enters  the  volume.  For  this  computation  we  return  to  the  notation 
used  in  section  2.1,  where  Xik  represents  the  components  of  stream  i  with  respect 
to  the  x,  y,  and  2;  axes,  indexed  by  variable  k  =  1,2,3,  respectively.  Similarly, 
let  lik  be  the  average  radiance  predicted  for  stream  i  entering  through  the  x,  y, 
and  2;  input  walls  of  the  volume,  again  indexed  by  variable  k.  And  let  Tik  be 
the  mean  transmittances  through  the  cell  for  energy  entering  through  the  A;th 
input  face  for  stream  i  (same  as  the  T  statistics  described  in  Eqs.  (43)  through 
(45),  except  that  the  ‘standard’  orientation  results  have  been  transformed  to 
correspond  to  the  correct  walls  for  the  stream  in  question). 

For  diffuse  stream  i  we  may  determine  the  total  flux  entering  the  cell  as: 

3 

E,n,i  ='^w,'Y^XikIik  ,  (100) 

k=l 

where  A  is  again  a  cubical  cell  edge  length,  and  where  we  weigh  the  energy 
according  to  the  solid  angle  associated  with  the  stream  {^Wi)  and  the  mean 
cosine  of  incidence  across  the  cell  wall  (xik)-  The  total  energy  transmitted 
across  this  cell  for  stream  i  will  similarly  be 

3 

Etran,i  —  ^  ^  ^iik  lik  Tik  A  .  (101) 

k=l 

If  we  then  assumed  a  conservative  scattering  condition,  we  would  expect  that 
the  scattered  energy  emitted  from  the  volume  element  would  account  for  the 
remaining  energy.  Call  this  amount  Ecs,i  for  the  conservative  scattered  energy. 
The  actual  energy  scattered  under  a  single  scattering  assumption  would  be 
w  Ecg^i,  where  ru  is  again  the  single  scattering  albedo.  Let 

M  3 

Ecs,z  =  R  ^  ^  Xjk  (1  -  Tjk)  =  C  k,  (102) 

j=l  k=l 

where  the  result  is  a  linear  function  of  R,  the  j  summation  is  over  all  output 
directions  from  the  cube,  R  is  the  appropriate  volume  averaged  radiance  for 
stream  i  that  produces  the  correct  amount  of  emitted  energy  in  the  conservative 
scattering  limit.  The  k  summation  is  over  the  output  faces,  and  thus  Xjk  is  the 
appropriate  cosine  of  the  jth  stream  with  the  A:th  output  face;  similarly,  Tjk  is 
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the  mean  transmittance  of  energy  exiting  the  A;th  face  via  the  jth  stream,  and 
jiji  is  the  cosine  of  the  angle  between  incident  (O^)  and  exiting  (Oj)  directions. 

Using  Eqs.  (100)  through  (102),  we  may  write 


(E 


in,i 


Efrariji 
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Z]  Xik  (1  —  Ttk)  Ilk 
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Ef 
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E{j^ji)  E  Xjk  (1 


k=l 


Tjk) 


(103) 


Using  this  value  for  the  volume  averaged  radiance  for  stream  z,  the  total  limiting 
path  radiance  for  all  streams  in  the  standard  discrete  ordinates  directions  O;  may 
then  be  determined  as. 


M 

k  =  w  (104) 

i  =  l 

However,  this  value  accounts  for  only  the  diffuse  scattering  component  of  the 
limiting  path  radiance.  To  this  value  we  must  add  the  black  body  component 
and  the  direct  component.  The  black  body  component  is  simply  the  B  term 
previously  described.  For  the  direct  component,  we  must  follow  a  similar 
procedure  to  that  described  for  the  diffuse  component,  but  the  derivation  must 
consider  the  delta  scattering  property  of  the  direct  radiance.  The  integration 
over  output  faces  for  the  conservative  scattering  calculation  involves  terms 
P(Oo,  where  Oq  is  again  the  unit  vector  indicating  the  direction  of 

propagation  of  the  direct  energy  component.  Second,  instead  of  describing  effects 
by  scattering  cell  class,  each  cell  must  be  individually  assessed  to  determine  the 
incident  flux  on  each  input  face  and  the  fractional  transmitted  energy  for  energy 
that  entered  through  each  face.  Lastly,  since  the  direct  energy  is  either  converted 
completely  into  diffuse  stream  energy,  transmitted  through  the  cell,  or  absorbed, 
there  is  no  additional  effect  of  direct  radiation  for  path  radiance  purposes. 

We  thus  have  a  complete  procedure  for  evaluating  the  limiting  path  radiance  for 
any  cell  in  any  direction,  and  thus  we  can  determine  the  apparent  brightness  of 
any  path  through  the  scattering  volume. 

4.4  Directional  Interpolation 

The  remaining  step  in  evaluation  of  the  limiting  path  radiance  is  a  directional 
interpolation.  The  reasoning  behind  this  step  is  that  instead  of  computing 
directly  based  on  Eq.  (104),  this  equation  requires  an  interpolation  over  all 
incident  directions  each  time  it  is  invoked.  An  alternative  method  of  evaluating 
would  be  to  translate  the  results  produced  at  the  discrete  ordinates  directions 
into  some  other,  higher  order  representation.  In  this  representation,  more 
efficient  methods,  using  a  smaller  set  of  nearest  neighbor  points,  could  be  used. 
One  such  possibility  is  to  generate  a  series  of  points  at  equally  spaced  meridians. 


81 


A  splined  fit  can  then  be  performed  along  each  meridian,  and  results  at  different 
meridians  can  then  be  further  splined.  Other  methods  are  possible  as  well;  for 
example,  one  may  desire  to  transform  from  the  set  of  discrete  ordinate  directions 
to  a  separate  set  of  new  directions  which  also  span  the  unit  sphere.  To  make 
this  transformation,  a  spherical  harmonics  method  is  invoked. 

Let  us  define  a  set  of  orthonormal  functions  on  the  unit  sphere  which  shall  be 
referred  to  as  even  and  odd  spherical  harmonic  functions: 


(t>)=  y  ^4^  ^  (2  ~  ^  cos(m(^),  (105) 


Of  1  (f  —  m  V 

ZUv.  Ct>)  =  2  Prifi)  Sin(m<^),  (106) 

where  /i  is  the  cosine  of  the  zenith  angle,  ^  is  an  azimuthal  angle,  f  >  m  >  0, 
and  where  is  the  associated  Legendre  polynomial,  which  can  be  given  in 

terms  of  Rodrigues’  formula  as  (Jackson  1975), 

(_^\m  Ji+m 

prw  =  ^4-  (1  -  -  !)'■  (107) 

These  functions  have  the  properties  of  orthogonality  when  integrated  over  dvr 
steradians: 

J  4)  =  (108) 

47r 

J  d9.Z1,^,{ji,  (f))Z^^(iJ.,  (f>)  =  -  S^o)]  (109) 

47r 

J dfiZf,^,{fi,  4)Z|^(/i,  4)  =  0,  (110) 

47r 

where  Sij  is  a  kronecker  delta  function  of  indices  i  and  j . 

These  functions  are  similar  to  those  used  by  Zardecki  (1995),  but  contain 
different  normalization  factors,  which  tend  to  somewhat  simplify  the 
mathematics.  Note  that  these  spherical  harmonics  functions  are  always  real 
valued.  This  is  significant  because  we  want  to  exclude  imaginary  results,  since 
we  are  dealing  with  strictly  real  valued  radiances. 

Because  of  these  properties,  a  function  of  direction  may  be  expanded  as: 

OO  i 

/(/i,  4)  =  E  E  n  (111) 

1=0  m=0 
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where  the  Ai^n  coefficients  are  determined  using, 
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im 


J  (p); 

47r 


(112) 


Aim  =  j  4i)I{ii,  i).  (113) 

47r 

However,  due  to  the  characteristics  of  the  discrete  ordinates  approach,  these 
coefficients  can  be  determined  using  the  summation  results: 


M 

=  (114) 

i=l 

M 

cbi)I,,  (115) 

i=l 

where  we  recall  from  section  2.2  that  M  (equals  N(N  +  2))  is  the  number  of 
streams  in  a  discrete  ordinates  model. 

If  one,  therefore,  desired  to  determine  the  radiance  results  (Ij)  in  a  series 
of  directions  Oj,  the  transform  mechanism  for  converting  from  the  sampled 
directions  li  to  the  output  directions  Ij  could  be  expressed  as  a  matrix 
multiplication: 
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is  obtained  by  rearranging  the  order  of  the  summations  in  Eqs.  (HI),  (114),  and 
(115),  and  Zp^^  =  (pi),  Z^^^  =  Z^°^(/i*,  (pi).  The  only  restriction  placed 

on  this  method  is  that  the  value  of  the  summation  limit  L  used  in  Eq.  (117) 
must  be  chosen  such  that  all  equations  of  form  (108)  through  (HO)  integrate 
properly  under  the  stream  expansion  chosen.  That  is,  we  require. 
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(120) 


47r  ^  \ 

^  ^  Wi  (f)i)  Z£^(jli,  (f)i)  ~  0. 

i=l 

This  requirement  leads  to  the  interesting  conclusion  that,  while  we  have  assumed 
that  if  M  =  iV  {N  +  2),  our  phase  function  expansion  was  valid  to  T  =  iV,  direct 
calculations  show  that  the  interpolation  criteria  above  are  only  valid  through 
L  =  N/2. 

Using  the  above  technique,  the  radiance  data  produced  by  the  radiative  transfer 
code  may  be  interpreted  for  any  desired  direction. 

In  summary,  we  may  conclude  that  transmittances  and  limiting  path  radiances 
for  individual  path  segments  can  be  computed  along  any  path  within  the 
modeled  scattering  volume,  leading  to  a  constructive  method  for  evaluating 
path  transmittances  and  integrated  path  radiances  for  arbitrary  geometries.  The 
remaining  aspects  of  propagating  information  for  image  rendering  purposes  are 
relating  the  image  geometry  to  a  given  path  within  the  scattering  volume  and 
relating  the  radiances  produced  by  the  scattering  model  to  pixel  data  placed  in 
the  image  file.  These  aspects  are  discussed  in  the  next  two  sections. 

4.5  Viewing  Geometries 

In  perspective  viewing,  the  origin  point  for  each  pixel  is  some  image  plane  within 
the  observer’s  optical  system.  The  simplest  such  system,  and  the  one  adopted 
here  for  our  case  study,  is  a  pinhole  camera  where  the  image  plane  is  placed 
perpendicular  to  the  primary  viewing  axis,  and  each  LOS  is  viewed  as  passing 
out  of  the  camera  through  the  center  of  the  lens  into  the  modeled  volume. 
In  orthographic  rendering,  the  origin  of  each  LOS  can  be  a  point  in  a  plane 
from  which  parallel  lines  of  sight  are  drawn  through  the  volume  of  interest.  In 
perspective  viewing,  the  lines  of  sight  naturally  diverge  after  passing  through 
the  ‘lens’  to  span  a  particular  solid  angle  view  of  the  volume  of  interest. 

Regardless  of  this  distinction,  each  LOS  will  trace  its  way  through  a  unique 
sample  of  the  modeled  scattering  volume  until  it  either  exits  the  volume  through 
the  top  of  the  modeled  region  (the  ‘roof’),  impacts  the  modeled  surface  (the 
‘floor’),  or  reaches  some  preselected  maximum  range  after  which  the  computation 
is  truncated.  This  latter  step  is  usually  necessary  to  control  computation  time 
in  an  application. 

Any  LOS  that  crosses  upward  through  the  roof  is  assumed  to  see  only 
background  atmosphere  with  a  brightness  characterized  by  a  directionally 
dependent  background  sky  color.  An  LOS  which  strikes  the  surface,  on  the  other 
hand,  can  be  represented  using  any  of  a  series  of  possible  models  to  include  a 
constant  reflectivity  surface,  user-designed  surface  objects,  reflectivity  mappings, 
etc. 
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4.6  Radiance  to  Pixel  Value  Conversions 

Once  means  of  representing  the  upper  and  lower  bounding  surfaces  have  been 
chosen,  it  will  be  necessary  to  convert  path  radiance  values  given  in  physical 
units  into  image  related  pixel  values.  This  process  involves  the  use  of  conversion 
factors. 

To  understand  where  AIM  derives  its  radiance  information,  we  begin  with 
another  AIM  model  input.  To  initialize  the  radiance  information  AIM  runs  an 
Electro-Optical  Systems  Atmospheric  Effects  Library  (Shirkey  et  al.  1987) 
(EOSAEL)  configured  version  of  the  Air  Eorce’s  MODTRAN  model  (Berk  et  al. 
1989)  to  determine  the  direct  (solar)  irradiance  and  directionally  dependent 
diffuse  radiances  of  the  sky  background  at  the  top  of  the  modeled  volume  in 
each  of  the  downward  welling  stream  directions.  Using  these  values  along  with 
the  interpolation  technique  described  in  section  4.4,  one  can  determine  the  sky 
radiance  in  any  direction.  Call  this  amount  Ru{^).  Due  to  the  nature  of 
MODTRAN,  Ru  will  only  be  a  function  of  direction,  not  position.  Once  an  LOS 
crosses  upward  out  of  the  volume,  using  the  integration  convention  described 
in  section  4.2,  we  would  simply  add  in  a  final  step  to  the  path  integration 
procedure  (call  it  iV  -|-  1),  where  s(y_|_^  =  oo,  =  0,  and  =  Ru-  Thus 

there  is  no  transmission  beyond  this  final  step,  and  Ru  becomes  the  background 
radiance.  Depending  on  how  high  above  the  surface  one  wanted  to  simulate 
embedded  (flying)  objects,  other  options  are  also  available  but  might  involve 
running  larger  simulation  volumes. 

Eor  downward  directed  LOS,  we  need  to  consider  lines  that  pass  through  the 
ground  plane.  To  accomplish  this  task,  we  need  to  know  the  properties  of  the 
surface  to  be  modeled  as  well  as  the  simulated  ambient  illumination  at  each 
surface  point.  To  evaluate  this  ambient  illumination,  consider  that  at  each 
surface  point  there  will  be  a  set  of  downward  directed  streams  which  can  be 
integrated  over  solid  angle  and  cosine  weighted  to  produce  a  net  diffuse  flux. 
The  direct  irradiance  can  be  computed  as  well,  resulting  in  a  net  irradiance  at 
the  base  of  the  modeled  volume  beneath  the  last  layer  of  modeled  cells.  Currently 
we  assume  the  surface  is  a  Lambertian  reflector.  Thus  the  reflected  radiance  can 
be  evaluated  by  dividing  the  net  incident  flux  by  vr  and  multiplying  by  a  surface 
reflectivity  a. 

The  BLITS  model  can  accept  positionally  varying  surface  reflectivity,  but  this  is 
not  the  general  usage  of  the  model.  Instead,  an  average  surface  reflectivity  is 
usually  applied  to  the  entire  surface.  Using  this  approach,  the  radiance  fields 
produced  can  be  translated  across  the  terrain,  reflecting  the  first  order  effects 
of  wind  motion  of  the  cloud  held.  Most  terrain  variations  will  average  out  over 
the  250  m  cell  lengths  involved  in  the  RT  calculations. 

Taking  this  approach,  assume  for  a  given  region  of  the  surface  that  there  is  a 
known  mean  reflectivity  d.  Assume  also  that  there  is  an  image  of  the  surface 
available.  Define  as  the  reflectivity  of  the  zero  pixel  value  and  um  as  the 
reflectivity  of  the  maximum  possible  pixel  value.  Also  assume  the  image  has  not 
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been  distorted  in  some  way,  such  as  by  a  gamma  correction.  Thus  the  image 
brightness  is  considered  to  vary  linearly  according  to  pixel  value,  and  that  this 
linearity  is  related  to  pixel  reflectivity. 

We  have  thus  discussed  means  of  representing  the  upper  and  lower  boundaries 
of  the  rendered  scene  using  previously  available  image  data  and  radiance 
information  from  MODTRAN.  The  following  section  completes  the  information  on 
rendering  techniques. 

4.7  Cloud  Rendering  Algorithm 

In  previous  sections  we  described  the  means  of  tracing  paths  through  the 
scattering  volume.  However,  as  is  often  the  case,  this  technique  is  not  without 
difficulties.  These  difficulties  arise  due  to  the  level  of  complexity  of  the  problem 
addressed.  For  a  tactically  significant  scenario,  we  would  like  to  treat  a  volume 
on  the  order  of  several  kilometers  in  each  horizontal  direction.  In  addition, 
we  need  to  extend  the  modeled  scattering  volume  from  the  surface,  where  the 
majority  of  tactically  significant  (from  an  Army  point  of  view)  lines  of  sight 
exist,  up  through  a  level  that  encompasses  the  majority  of  significant  cloud 
features.  This  extension  vertically  is  necessary  because  cloud  shadowing  can 
have  significant  effects  on  the  illumination  of  the  target  and  also  the  LOS, 
significantly  affecting  path  radiance. 

The  difficulty  with  extending  the  scenario  to  encompass  a  tactically  significant 
volume  is  that  once  it  has  been  sufficiently  extended,  the  small  feature  data  in  the 
clouds  themselves  often  cannot  be  adequately  characterized  during  the  radiative 
transfer  calculation.  In  our  typical  scenario  a  volume  with  8  km  horizontal  by 
4  km  vertical  dimensions  is  modeled  using  A  =1/4  km  per  cell.  But,  this  cell 
size  is  often  the  size  of  actual  clouds.  Thus,  a  cloud  may  be  modeled  by  only  a 
single  cell.  This  can  lead  to  significant  problems  when  attempting  to  visualize 
such  a  cloud.  Without  additional  information  the  cloud  will  appear  as  a  cube. 

To  avoid  this  problem,  an  observation  made  by  Hoock  and  Giever  (1989)  was 
utilized.  This  observation  was  that  limiting  path  radiance  varies  more  slowly 
than  path  radiance,  since  path  radiance  depends  on  the  specific  density  structure 
of  the  medium,  but  limiting  path  radiance  reflects  the  presence  of  available 
illumination.  In  most  clouds,  the  illumination  is  due  to  scattered  energy,  which 
becomes  rather  uniform  due  to  the  highly  scattering  nature  of  most  cloud 
aerosols.  This  suggests  that  the  radiance  calculations  in  the  radiative  transfer 
model  can  be  run  at  a  coarser  resolution  than  is  needed  by  the  density  map  of 
the  medium. 

To  generate  the  3D  extinction  mapping  of  the  scattering  medium,  the  AIM 
preprocessing  stage  runs  the  Air  Force’s  Cloud  Scene  Simulation  Model  (CSSM). 
At  its  heart  this  model  contains  a  fractal  density  generator  and  fractal  statistics 
for  various  cloud  types.  From  the  AIM  input  interface,  a  user  can  select 
characteristics  of  the  cloud  held  to  be  generated,  including  the  base  heights. 


layer  thicknesses,  cloud  types,  and  cloud  fractions  of  up  to  three  layers  of  non¬ 
cumulus  clouds  and  one  cumulus  cloud  layer.  AIM  then  runs  CSSM  at  four  times 
the  resolution  of  that  required  by  the  BLITS  RT  code.  The  full  resolution  data 
generated  by  CSSM  are  then  transformed  to  extinction  information  which  can  be 
used  in  visualizations.  The  data  are  then  also  coarsened  using  a  spatial  averaging 
function  for  use  by  BLITS. 

A  rendering  code  was  then  developed  to  access  the  output  from  the  RT 
code  and  the  3D  extinction  mapping.  The  various  techniques  for  generating 
pixel  modifications,  angular  interpolation,  and  point-to-point  integrations 
were  employed.  The  model  developed  was  tailored  to  generate  orthographic 
renderings  of  the  scattering  volume.  The  model  entails  two  separate  calculation 
stages.  In  the  first,  a  series  of  lines  of  sight  are  traced  through  the  volume.  In 
the  second,  these  results  are  used  as  interpolation  points  for  generating  values 
at  every  pixel  position. 

In  the  first  phase,  the  1/16  km  resolution  CSSM  output  is  used  in  combination 
with  the  1/4  km  resolution  BLITS  model  output.  For  a  given  observer  location 
on  the  outside  of  the  modeled  volume,  there  will  be  at  most  three  outside  walls  of 
the  volume  that  are  visible.  Normally,  these  walls  cover  a  simulated  8  kmx8  km 
footprint  on  the  ground  and  are  4  km  high.  They  are  sampled  at  the  same 
1/16  km  resolution  as  output  from  the  CSSM  code.  Thus,  on  a  vertical  side  wall 
8  km  wide  and  4  km  high,  there  are  128x64  sample  lines  of  sight  traced  through 
the  volume.  On  a  top  or  bottom  wall  there  will  be  128x128  sample  points.  Each 
LOS  is  traced  to  determine  net  transmission  and  path  radiance  for  each  of  three 
color  channels. 

Following  these  LOS  calculations,  the  program  uses  this  data  set  to  evaluate  the 
results  for  each  pixel.  A  given  pixel  will  be  characterized  by  a  path  through  the 
volume.  A  Gaussian  weighting  scheme  is  used  to  interpolate  the  results  for  pixels 
falling  between  LOS  calculations  made  in  the  previous  phase.  This  technique 
provides  a  smooth  transition  between  effects  calculated  for  the  sample  points, 
as  characterized  by  the  width  of  the  Gaussian  blending  functions.  Gare  must 
be  taken  in  selecting  an  appropriate  width  parameter,  since  a  width  parameter 
that  is  too  small  will  produce  grainy  results  where  only  the  closest  sample  LOS 
calculation  can  influence  the  result.  A  width  parameter  which  is  too  large  will 
cause  clouds  to  appear  too  smooth  and  can  greatly  slow  down  the  interpolation 
algorithm. 

For  a  rendering  algorithm,  one  must  also  describe  the  appearance  of  the  surface 
or  sky  backgrounds.  We  have  created  a  separate  file  that  lists  the  incident 
direct  and  diffuse  radiances  present  at  each  point  beneath  the  scattering  volume 
for  each  sensor  color  channel.  Let  us  call  the  combined  channel  irradiance 
RljA’  ''^here  i  and  j  represent  the  horizontal  position  across  the  terrain,  and 
k  represents  the  color  channel  index.  As  an  example,  a  color  sensor  receiving 
data  on  red,  green,  and  blue  channels  has  k  values  of  1,  2,  and  3.  Let  us  also 
define  surface  color  channel  reflectivities  resulting  in  reflected  energy, 

Oii,j,kVij^k  at  each  (z,  j)  point  along  the  surface  and  for  each  color  channel  k. 
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For  upward  directed  paths,  there  are  corresponding  values  to  be  applied  for  sky 
radiance  color  channel  based  on  the  description  given  in  section  4.6.  The  last 
step  is  then  to  apply  the  transmittance  and  path  radiance  effects  interpolated 
from  the  initial  data  set  to  the  background  pixel  values  just  determined. 

The  results  of  these  processing  stages  are  a  full  set  of  real  number  results  at  each 
pixel  position  for  the  color  channel  radiance  results.  The  final  step  necessary  is 
to  rescale  these  results  such  that  they  are  remapped  into  a  range  that  does  not 
exceed  the  maximum  pixel  value  possible  in  any  color  channel.  To  accomplish 
this  step  requires  two  passes  through  the  pixel-based  data  set.  The  first  step 
involves  determining  the  maximum  pixel  value  at  any  channel  and  then  using 
this  as  a  normalizing  factor  in  the  second  stage. 

Sample  outputs  from  this  algorithm  are  seen  in  figures  44  through  46.  The 
scenario  modeled  uses  a  green,  20-percent  reflective,  terrain  surface  for  the  base 
of  the  volume.  The  cloud  held  is  that  for  a  35  percent  cumulus  cloud  cover.  The 
cloud  held  has  been  computed  such  that  it  is  horizontally  periodic.  Thus  clouds 
that  begin  on  one  side  or  corner  of  the  volume  wrap  around  to  the  opposite 
edge.  This  technique  was  deemed  preferable  to  either  (a)  modeling  very  large 
volumes  or  (b)  modeling  clouds  which  simply  cut  off  at  the  boundaries  of  the 
volume. 

Orthographic  views  of  the  scattering  volume  are  shown  in  figures  44  through  46 
as  seen  from  the  top,  side,  and  below.  An  afternoon  time  is  modeled,  resulting 
in  illumination  on  the  western  side  of  the  clouds.  Dark  regions  appear  on  some 
clouds  because  the  horizontal  periodicity  results  in  some  clouds  on  the  eastern 
edge  shading  clouds  near  the  western  edge.  In  particular,  note  the  shadowed 
cloud  tops  in  the  center  of  figure  44.  Figure  45  reveals  that  the  cloud  tops  of  the 
clouds  on  the  eastern  edge  (in  the  background  in  figure  45)  are  actually  taller 
than  the  clouds  in  the  center.  They  thus  shadow  the  clouds  in  the  center  due 
to  the  horizontal  periodicity.  Cloud  shadows  also  result  in  surface  brightness 
variations,  which  show  up  better  on  the  computer  screen  than  on  hard  copy. 

4.8  Fast  VIEW  Data  Formats 

The  AIM  codes  utilize  a  data  formatting  and  compression  technique  developed, 
but  not  formally  documented,  in  1994  for  the  Defense  Modeling  and  Simulation 
Office’s  (DMSO)  Environmental  Effects  for  Distributed  Interactive  Simulations 
(E^DIS).  In  that  effort,  an  attempt  was  made  to  generate  common,  useable,  yet 
compact  data  sets.  The  name  established  for  these  data  sets  was  East  VIEW, 
reflecting  the  relationship  between  their  source,  the  Weather  and  Visualization 
Effects  for  Simulations’  (WAVES)  model  VIEW  and  an  attempt  to  utilize  data 
output  by  that  model  in  a  rapid  manner.  We  now  provide  the  documentation 
formally  for  that  format  in  this  section  in  hopes  that  others  can  design  uses  that 
provide  faster  rendering  of  images  using  these  sets.  Originally,  these  formats 
consisted  of  only  two  data  set  types.  We  have  expanded  these  to  three  here 


Figure  44.  Visualization  of  a  35  percent  cloud  cover  cumulus  layer.  View  is  from 
directly  above  using  an  orthographic  rendering  technique.  North  is  at  the  top  of  the 
figure. 


Figure  45.  Visualization  of  a  35  percent  cloud  cover  cumulus  layer.  View  is  from  the 
western  side  of  the  volume  looking  east  using  an  orthographic  rendering  technique. 
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Figure  46.  Visualization  of  a  35  percent  cloud  cover  cumulus  layer.  View  is  from 
directly  below  using  an  orthographic  rendering  technique.  North  is  at  the  top  of  the 
figure. 

to  reflect  the  need  for  a  higher  density  extinction  map  useable  by  a  rendering  engine. 

These  data  types  now  consist  of  the  file  types:  EX.OUT,  IL.OUT,  and  FV.CPD. 
The  first  provides  information  on  the  density  map  of  the  volume.  The  second  is  a 
surface  illumination  map.  The  third  provides  volumetric,  directional,  and  spectral 
radiance  information. 

4.8.1  EX.OUT 

The  extinction  file,  EX.OUT,  provides  a  3D  spectral  mapping  of  the  scattering 
volume's  extinction  coefficient  mapping.  As  explained  in  the  previous  section,  the 
RT  code  is  run  on  a  coarse  grid  to  conserve  memory  and  time.  This  information  is 
then  coupled  to  a  higher  resolution  map  of  the  extinction.  The  EX.OUT  file 
provides  that  map. 

The  EX.OUT  file  itself  consists  of  two  sections,  a  header  containing  information 
about  the  volume  and  spectral  resolutions,  and  the  data  itself.  A  typical  header 
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is  seen  below: 


EXTC 

1.0000 

XDIS 

0.0625 

MXX 

128.0000 

YDIS 

0.0625 

MYY 

128.0000 

ZDIS 

0.0625 

MZZ 

64.0000 

MCC 

1.0000 

These  lines  represent  the  following  information:  The  first  line  merely  identifies 
the  file  type  as  an  EX .  OUT  file;  the  second,  fourth,  and  sixth  lines  indicate  the  cell 
resolutions  in  kilometers;  the  third,  fifth,  and  seventh  lines  indicate  the  number 
of  cells  in  the  three  directions;  and  the  last  line  indicates  the  number  of  spectral 
channels.  The  data  body  which  follows  the  header  section  is  listed  according 
to  spectral  channel,  x,  y,  and  z  in  order  from  most  to  least  rapidly  varying  cell 
index.  In  our  convention,  x  is  always  increasing  toward  the  east,  y  increases 
toward  north,  and  z  is  a  vertical  variable  measured  as  height  above  the  surface. 

4.8.2  IL.OUT 

The  surface  illumination  file  type,  IL  .  OUT,  is  similar  to  the  extinction  file  except 
that  there  are  only  two  dimensions  in  the  data  and  there  is  more  than  one  data 
element  per  line  in  the  body  of  the  output.  Here  is  a  listing  of  a  typical  header: 


ILUM 

1.0000 

XDIS 

0.2500 

MXX 

32.0000 

YDIS 

0.2500 

MYY 

32.0000 

MCC 

3.0000 

DZEN 

49.7397 

DAZI 

173.8699 

The  data  listed  for  this  file  are  the  direct  and  diffuse  fluxes  for  each  surface 
region.  The  data  are  organized  in  the  same  format  as  the  EX .  OUT  results: 
spectral  channel,  x,  and  y  in  order  of  nesting  from  deepest  to  outer  (y).  Note  that 
these  results  are  produced  in  some  cases  by  weighing  the  computed  radiances  at 
the  surface  according  to  their  incident  cosine,  solid  angle  for  diffuse  quantities, 
and  spectral  weighting  function. 

As  an  example  of  a  spectral  weighting  function,  consider  the  human  eye. 
Middleton  (1952)  discusses  the  characteristics  of  the  eye’s  photopic  and  color 
spectral  response  curves.  Let  these  be  characterized  by  functions  Sc{k)^  where  k 
is  the  wavenumber  and  c  is  the  color  channel  index.  These  curves  are  normalized 
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such  that  Sc(k)  dk  =  1.  Hence,  an  input  signal  of  E(x,  y,  k)  in  either  diffuse 
or  direct  flux  produces 


To  tailor  these  responses  for  use  with  image  rendering  products,  we  typically 
convert  the  Ic  values  to  NTSC  RGB  metrics  using  the  transform  i?  =  iV/, 
where  N  is  the  matrix. 


N 


1.910  -0.533  -0.288 

-0.985  2.000  -0.028 

0.058  -0.118  0.896 


(122) 


Users  should  note  that  in  this  form  the  effects  of  sensor  weighting  have  already 
been  accounted  for  once.  In  using  these  results  for  color  scene  rendering,  one 
needs  to  use  terrain  reflectance  coefficients  that  do  not  themselves  also  contain 
sensor  weighting  factors  or  the  sensor  effects  would  be  applied  twice. 

4.8.3  FV.CPD 

The  final  output  form  is  the  Fast  VIEW  compressed  file  type.  This  file  consists 
of  three  main  sections:  a  header,  a  series  of  directional  templates,  and  the  main 
data  listing.  Of  these,  the  first  section  is  virtually  identical  to  the  EX .  OUT  style 
header.  An  example  is  given  below: 


FSVU 

1.0000 

XDIS 

0.2500 

NNXX 

32.0000 

YDIS 

0.2500 

NNYY 

32.0000 

ZDIS 

0.2500 

NNZZ 

16.0000 

STDR 

0.0625 

NNCC 

3.0000 

NDCL 

128.0000 

Note  that  the  differences  include  the  file  type  identifier  (FSVU  in  this  case),  the 
STDR  identifier,  and  the  NDCL  identifier.  The  STDR  identifier  defines  the  standard 
range  used  in  the  computations  of  the  transmission  factors  contained  in  the 
general  output  section.  The  NDCL  identifier  defines  the  number  of  directional 
templates  included  in  the  template  section. 

Each  template  currently  consists  of  43  items  of  data:  the  template  index  number 
followed  by  42  directional  values.  These  values  range  from  zero  to  32767  (2^^  —  1) 
and  describe  the  directionally  varying  nature  of  the  limiting  path  radiance.  To 
use  these  results  an  accessing  algorithm  has  been  written  in  CH — In  this 
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algorithm,  the  set  of  templates  is  described  by  a  matrix  Dij^  where  i  is  the 
template  index  and  j  is  the  direction.  In  the  main  positional  data  listing,  for 
each  location  and  sensor  spectral  channel,  a  single  number  is  used  to  indicate 
the  template  to  use  to  describe  limiting  path  radiance  directional  variability, 
z(a;,  y,  2;,  c).  Also,  a  single  number  is  used  to  define  the  maximum  limiting 
path  radiance  in  any  direction  at  that  point,  Lmix-,y-,  z^c).  To  determine  the 
actual  limiting  path  radiance  in  one  of  the  canonical  directions  j,  one  evaluates 
Lm(x,y ,  z,  c)  D[i(x,y ,  z,  c),  j]/327Q7 .  Intermediate  directions  can  be  evaluated 
by  any  of  a  number  of  interpolation  techniques.  The  CH — h  routine  mentioned 
has  the  full  description  of  each  of  the  42  directions  and  means  of  interpolation. 

Following  the  template  listing  is  the  main  section  of  the  data.  This  listing 
consists  of  data  in  a  different  order  from  the  EX .  OUT  listing:  x,  y,  z,  color  channel, 
in  order  from  outer  to  inner  loops.  Each  data  line  contains  four  elements:  the 
spectral  extinction  coefficient,  a  transmittance  factor,  a  maximum  limiting  path 
radiance,  and  a  directional  template  index  number.  The  transmittance  factor  is 
generated  based  on  the  path  length  value  (in  km)  given  in  the  header,  which  can 
be  used  to  step  through  the  cloud  field.  In  the  current  implementation,  STDR  is 
always  set  at  0.0625  km. 
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5.  Discussion  and  Conclusions 


In  the  previous  sections  we  have  shown  a  method  for  producing  radiometrically 
correct  three-dimensionally  inhomogenous  cloud  scenes  using  a  cell-surface- 
based  DOMRT  code.  This  code  has  the  advantage  that  energy  conservation  is 
imposed  automatically.  Also,  interpretation  of  diffuse  radiance  and  scattering  in 
any  direction  is  obtained  using  spherical  harmonics  expansions  and  the  discrete 
ordinates  choices  of  stream  directions.  Diffuse,  direct,  and  blackbody-emitted 
energies  are  all  treated  in  a  consistent  manner.  This  method  does  not  share 
the  disadvantages  of  Fourier  methods  regarding  memory,  nor  does  it  share  the 
limitations  of  the  original  DOM  in  regard  to  treating  only  low-density  media, 
nor  the  diffusion  approximation’s  limitation  to  dense  and  isotropically  scattering 
media.  Finally,  physically  realistic  cloud  scene  rendering  is  possible  using  results 
determined  by  the  RT  routine. 

The  performance  of  these  models  using  current  computer  resources  is  hopeful. 
Using  an  R4400  chip  SGI,  the  runtime  for  the  radiative  transfer  routine 
spanning  seven  wavebands  across  the  visible  spectrum  is  on  the  order  of  1 
hour.  Postprocessing  of  the  outputs  to  compress  the  information  to  color 
channels  and  allow  for  rapid  directional  interpolation  requires  about  1  additional 
hour.  Rendering  of  images  currently  takes  approximately  3  to  5  min  per  frame. 
However,  considering  that  no  parallelization  has  been  attempted  on  these  codes 
and  much  faster  machines  are  becoming  available,  it  is  possible  that  true  3D 
color  simulations  of  radiometrically  correct  cloud  fields  are  possible  in  the  near 
future.  The  compression  techniques  employed  are  also  useful  in  that  they  tend 
to  limit  the  amount  of  memory  that  must  be  dedicated  to  running  such  LOS 
calculations. 

This  method  has  been  integrated  with  an  intuitive  scenario  generator  within 
AIM.  This  section  of  the  code  allows  the  characterizing  of  a  particular 
scenario  using  a  minimum  of  parameters  while  automatically  processing 
all  the  needed  information  to  allow  the  rendering  of  the  scenario  volume 
by  the  rendering  engine.  Intermediate  output  files  needed  to  render 
the  volume  from  different  directions  are  archived  to  facilitate  additional 
simulation  tasks.  We  have  implemented  a  version  of  this  software  within  the 
DoD  DMSO-sponsored  Master  Environmental  Library  program  (http://www- 
mel.nrlmry.navy.mil/),  permitting  the  on-line  access  and  delivery  of  such 
files.  We  have  also  established  a  web  site  providing  color  examples 
of  the  images  that  can  be  created  using  this  technique  (http://www- 
mel . arl .mil/MEL/rad/rO .html). 
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AIM  allows  a  user  to  produce  any  number  of  different  cloud  scenarios  for 
visualization  purposes  and  provides  choices  in  terms  of  cloud  representation 
within  the  RT  code.  This  code  accesses  databases  of  MODTRAN  horizontally 
homogeneous  background  aerosols  profiles.  The  EOSAEL  configured  version 
of  MODTRAN  is  itself  run  several  times  to  initialize  the  volume  upper  boundary 
with  information  regarding  downwelling  diffuse  stream  radiance  and  direct 
illumination  as  well  as  being  used  to  compute  vertical  components  of  background 
molecular  absorption  coefficients.  It  also  calls  the  Air  Eorce  CSSM  code  to 
generate  a  cloud  density  field. 

AIM  uses  a  solar/lunar  calculation  routine  to  determine  whether  the  Sun  or  the 
Moon  is  the  source  of  direct  illumination.  Whenever  the  Sun  is  above  the  horizon 
it  is  used  as  the  direct  source.  If  the  Sun  is  below  the  horizon  and  the  Moon 
is  above  the  horizon,  then  the  Moon  is  used  as  the  direct  source.  Appropriate 
changes  to  inputs  parameters  used  in  the  calls  to  the  MODTRAN  routine  are  used 
to  reflect  the  conditions  to  be  simulated. 

AIM  also  uses  information  gleaned  from  the  PFNDAT  series  of  aerosols  contained 
in  EOSAEL.  This  series  of  aerosols  covers  the  majority  of  surface-based 
precipitation  types  (not  treated  under  the  current  version  of  AIM),  hazes,  fogs, 
and  munition  smokes.  In  addition  to  these  aerosols,  a  sequence  of  upper- 
level  cloud  aerosols  were  run  through  the  AGAUS  Mie  scattering  phase  function 
calculation  routine  to  produce  PFNDAT-clone  outputs.  The  PFNDAT  files  and 
the  clone  results  were  then  processed  through  a  routine  to  analyze  the  optimal 
Legendre  polynomial  representation  for  the  24-stream  model  using  the  technique 
described  in  the  appendix. 

The  issue  of  cloud  and  haze  visualization  is  a  held  of  intense  current  interest. 
Though  many  current  and  planned  simulations  attempt  to  treat  cloud  effects, 
these  seldom  treat  clouds  in  radiometrically  correct  fashion  or  consider  fully 
three-dimensional  influences.  More  often  than  not,  the  appearance  of  clouds 
are  handled  as  merely  pictures  that  have  been  artificially  embedded  in  a  scene, 
without  regard  to  their  feedback  interactions  with  scene  features  such  as  via 
cloud  shadow  effects  or  via  limiting  path  radiance  variations  produced  by  shafts 
of  sunlight  extending  through  gaps  in  the  cloud  layers. 

The  codes  described  here  are  far  more  robust.  In  the  future,  the  processing  speed 
of  these  codes  can  be  greatly  enhanced  using  parallel  processing  techniques,  novel 
rendering  techniques,  and  increased  computer  memory. 
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Acronyms  and  Abbreviations 


ID 

3D 

AIM 

ASL 

AFGL 

ARL 

AGAUS 

AGSNOW 

BLIRB 

BLITS 

BRDF 

GSSM 

DISORT 

DOM 

DOMRT 

DMSO 

DoD 

E^DIS 

EOSAEL 

EOSAEL92 

EE 

HOPE 

IHMG 

IID 

LLS 

LOS 

M&G 

MODTRAN 

RMS 

RMSE 

RT 

RVs 

SAA 

WAVES 


One- dimensional 

Three-dimensional 

Atmospheric  Illumination  Module 

U.S.  Army  Atmospheric  Sciences  Laboratory 

U.S.  Air  Eorce  Geophysics  Laboratory 

U.S.  Army  Research  Laboratory 

August  Miller’s  Mie  Scattering  Gode 

Special  version  of  AGAUS  for  treating  SNOW  cases 

Boundary  Layer  Illumination  and  Radiance  Balance  model 

Boundary  Layer  Illumination  and  Transmission  Simulation 

Bi-directional  Reflectance  Distribution  Eunction 

Air  Eorce’s  Gloud  Scene  Simulation  Model 

Discrete  Ordinates  Radiative  Transfer  routine 

Discrete  Ordinate  Method 

Discrete  Ordinate  Method  approach  to  Radiative  Transfer 
Defense  Modeling  and  Simulation  Office 
Department  of  Defense 

Environmental  Effects  for  Distributed  Interactive  Simulations 

Electro-Optical  Systems  Atmospheric  Effects  Library 

1992  Release  of  EOSAEL 

Einite  Element  Technique  within  DOM 

Henyey-Greenstein  Phase  Eunction 

In-House  Monte  Garlo 

Independent  Identically  Distributed 

Log  Least  Squares 

Line  of  Sight 

McKee  and  Gox 

Air  Eorce  Moderate  Resolution  RT  Model 

Root  Mean  Square 

Root  Mean  Square  Error 

Radiative  Transfer 

Random  Variables 

Small  Angle  Approximation 

Weather  and  Atmospheric  Visualization  Effects  for  Simulations 
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Appendix  A 

The  Scaling  Transformation 
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In  this  appendix,  we  describe  the  approach  used  to  make  a  scaling  transformation 
on  the  form  of  the  phase  function  used  in  the  radiative  transfer  (RT)  code. 
This  transform  modifies  the  coefficients  used  to  describe  the  extinction,  single 
scattering  albedo,  and  Legendre  polynomial  expansion  coefficients  of  a  phase 
function,  in  the  process  removing  the  effects  of  the  aerosol  forward  scattering 
peak  from  the  radiative  transfer  equation. 

These  modifications  are  necessary  because  phase  functions  often  have  forward 
peak  values  on  the  order  of  hundreds  or  thousands.  It  is  impossible  to  adequately 
simulate  these  peaks  with  low  order  Legendre  expansions.  Thus,  it  is  a  common 
practice  to  model  the  forward  peak  region  with  a  S  phase  function.  By 
introducing  this  approximation,  some  of  the  energy  is  treated  as  unscattered. 

The  method  developed  here  consists  of  a  log-least-squares  (LLS)  approach  to 
evaluate  the  forward  scatter  fraction  (/)  and  an  empirical  model  for  modifying 
the  Legendre  expansion  coefficients,  c^.  The  LLS  approach  models  the  phase 
function’s  forward  peak  using  a  Gaussian  phase  function  adding  to  the  low 
order  Legendre  expansion  used  to  characterize  the  remainder  of  the  angular 
dependence. 

In  developing  the  empirical  form  for  the  c^,  we  consider  several  forms  for  the 
least  squares  approach.  It  is  shown  that  the  (5-M  solution  for  the  c^,  as  a  function 
of  /  and  the  Legendre  expansion  coefficients  of  the  phase  function  derives,  from 
a  standard  least  squares  analysis.  We  then  show  the  similarity  between  actual 
phase  functions  and  the  Henyey-Greenstein  single-parameter  phase  function. 
This  similarity  is  used  to  argue  that  a  further  correction  is  needed  to  the  ci 
equation,  leading  to  the  empirical  form.  The  results  of  using  this  technique 
have  previously  been  shown  in  figures  17  through  20  in  chapter  3,  which  reveal 
that  we  are  able  to  outperform  the  (5-M  technique  in  most  cases. 

A.l  Background 

As  described  in  the  main  contents  of  this  report,  the  discrete-ordinates  approach 
to  solving  the  radiative  transfer  equation  is  a  valuable  tool  in  studying 
scattering  in  three-dimensional  media.  A  main  feature  of  this  technique 
involves  representing  phase  function  scattering  in  a  way  that  is  compatible 
with  the  discrete  ordinates  methods.  To  accomplish  this  representation  we 
first  produces  estimates  of  the  scattering  phase  functions  of  various  natural 
aerosols.  Tofsted  et  al.  (1997)  generated  a  series  of  such  functions  for  near 
surface  haze  and  precipitation  aerosols  based  on  the  Mie  scattering  algorithm 
developed  by  Miller  (1983),  that  incorporated  the  improved  forward  peak 
handling  technique  of  Lentz  (1976).  Upper  air  cloud  aerosols  were  also  treated 
using  the  Khrgian-Mazin  particle  size  distribution  parameterizations  adopted  by 
Low  and  O’Brien  (1987). 
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These  phase  functions  must  then  be  represented  as  Legendre  expansions,  and 
often  for  greater  accuracy  these  expansions  are  modified  according  to  the 
delta  correction,  also  known  as  a  scale  transformation  (Potter  (1970),  Joseph 
et  ah  (1976),  Wiscombe  (1977),  Schaller  (1979),  McKellar  and  Box  (1981), 
LeNoble  (1985)).  Techniques  developed  under  these  previous  applications 
normally  only  considered  ID  RT  codes  where  large  order  stream  models  could 
be  instituted.  In  these  applications  little  consideration  was  given  to  ensuring 
that  the  phase  function  remained  nonnegative  at  all  scattering  angles  (with  the 
exception  of  Fiveland  (1984)  when  he  limited  the  asymmetry  parameter  g  to 
the  range  0<(/<l/3in  his  two-term  expansion)  because  in  most  situations 
the  scattering  model  order  was  high  enough  that  it  would  automatically  be 
nonnegative.  Except  for  Fiveland,  previous  studies  never  even  considered  that 
the  order  of  the  stream  model  would  be  too  low.  But  since  our  abilities  are 
already  being  taxed  due  to  CPU  time,  memory,  spectral,  spatial  resolution,  and 
spatial  domain  concerns,  it  only  makes  sense  that  we  should  have  to  implement 
our  model  under  particularly  stressful  conditions.  We  therefore  require  a  model 
which  is  more  robust  at  low  order  expansions  than  has  been  available  heretofore. 

To  avoid  the  problems  of  negative  radiance  calculations,  the  traditional  delta 
correction  approach  has  been  revised.  In  the  past,  various  empirical  arguments 
were  made  for  determining  the  value  of  the  parameter  /.  Joseph  et  al.  (1976) 
used  f  =  .  Liou  (1992)  proposed  using  a  constant  /  for  a  given  scattering 

order.  McKellar  and  Box  (1981)  recommended  various  functions  of  the 
Legendre  coefficients.  Most  of  these  recommendations  depend  on  the  number 
of  coefficients  selected.  Wiscombe  (1977)  set  /  to  the  next  Legendre  expansion 
coefficient  beyond  the  range  of  coefficients  modified  in  the  expansion.  But,  in 
each  of  these  cases,  /  is  defined  based  on  an  empirical  rule.  In  what  follows  we 
base  /  on  a  LLS  analysis,  but  determine  the  ci  expansion  coefficients  based  on 
an  empirical  rule. 

To  help  define  the  degree  of  success  possible  with  this  method  we  have  selected 
one  of  the  Shettle  and  Fenn  (1979)  advection  fog  aerosols  as  our  sample  case  to 
model.  This  aerosol  has  a  modified  gamma  particle  size  distribution  with  mode 
radius  of  ro  =10  gm  and  a  and  7  parameters  of  1  and  6,  respectively.  The 
scattering  properties  of  this  phase  function  at  different  wavelengths  are  modeled 
using  a  Mie  scattering  code. 

A. 2  Phase  Function  Representation 

As  in  the  main  section,  the  phase  function  P{g)  is  represented  as  a  function  of 
the  cosine  of  the  scattering  angle  (g)  only.  This  presupposes  that  the  aerosol 
droplets  are  either  spherical  in  shape  or  exhibit  random  orientation.  P{g)  is 
normalized  according  to 


2n 


d(f)  /  P[cos(0)]  sin(6)d6  =  2tt  j  P(g)dg  =  1. 


-1 


(A.l) 
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This  phase  function  can  be  exactly  represented  by  an  infinite  Legendre  series, 


P{li)  =  Xi  =  2Tr  I  at  =  ^  ,  (A. 2) 


where  Pi{fJ,)  is  the  Legendre  polynomial  of  order  £.  Unfortunately  we  do  not  have 
an  infinite  number  of  streams  by  which  to  characterize  the  scattering  effects,  so 
we  must  find  some  suitable  way  to  truncate  this  expansion. 

Under  the  delta  correction  approach  (Joseph  et  ah  (1976),  Wiscombe  (1977)), 
the  phase  function  is  normally  approximated  by  a  (5  function  in  the  forward 
direction  and  a  Legendre  expansion  of  order  L, 

Pe{p)  =  f  P  {I  -  f)'Y^aiciPi{ji).  (A. 3) 


Here  the  8  is  divided  by  2tt  after  van  de  Hulst  (1980b),  which  assumes  the  full 
peak  falls  within  the  integral  bounds. 

Under  the  classic  (5-M  approach  (Wiscombe  1977),  one  assumes  /  =  Xi_|_i.  Since 
the  Legendre  expansion  of  8{jj.  —  l)/(27r)  has  coefficients  Xi  =  1,  for  all  7,  if  one 
equates  the  first  L  +  2  terms  of  the  Legendre  expansion  of  P(/i),  in  Eq.  (A. 2), 
to  those  of  Pe{p),  in  Eq.  (A. 3),  one  obtains  the  equations. 


whence. 


a^Xi  =  at  f  X  I  +  ai{l  -  f)ci,  7  =  0, ...,  L  +  1, 


C£ 


Xi-f 
1-/  ’ 


7  =  0,...,L  +  1. 


(.4.4) 


(.4.5) 


Of  course,  since  Xq  =  1  by  definition,  we  will  always  have  cq  =  1.  And,  as 
noted,  since  under  the  (5-M  approach  /  =  Xi_|_i,  Ci_|_i  =  0.  This  expansion  is 
thus  purported  to  apply  to  an  order  L  +  1  expansion,  although  one  only  retains 
terms  out  to  L. 


Using  the  model  of  Pe{p)  above,  we  can  determine  its  effect  on  the  overall 
equation  of  transfer.  Let  us  rewrite  the  equation  of  transfer  using  the  phase 
function  estimate,  as 


O  •  V/(U,  O)  +  (7  /(U,  O)  =  aw  /  /(U,  O^)  Pe{^  ■  +  (7  B(f,  A).  (A. 6) 

J  47r 

Due  to  the  8  function  in  Pe{p),  the  original  scattering  source  term  can  be  written 
as. 


J{Q,r)  =  fawI{f,Q)  +  {l-f)aw  /(r.  O')  P'(0  •  O')  O',  (A.7) 

J  i-K 
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where  is  the  modified  Legendre  expansion: 


L 

P'ip)  =  Y.  aiciPi{ii).  (A. 8) 

The  RT  equation  can  then  be  rewritten  as, 

O  •  V/(r,  O)  +  a'  I(f,  O)  =  a'w'  /  /(r,  O')  P'(0  •  O')  O'  +  a'  B' ,  (A.9) 

J  47r 

which  has  the  same  form  as  the  original  equation  of  transfer,  except  that  we 
have  introduced  the  variables 


.7'  =  <t(1  -  /ro). 


(,4.10) 


and 


,  _  {I-  f)w 
B'  =  {1  —  w')b, 


(A.ll) 

(A.12) 


where  we  note  that  a'(l  —  w')  =  a{l  —  zu)  such  that  the  emittance  term  is 
unaffected  by  this  transformation. 

This  new  equation  simulates  the  forward  scattered  radiation  not  by  directly 
calculating  the  results,  but  rather  by  treating  forward  scattered  energy  as 
unscattered.  That  is,  the  angular  resolution  achieved  by  using  a  smaller  number 
of  streams  is  insufficient  to  resolve  the  difference  between  forward  scattered 
and  unscattered  radiation.  The  resulting  transfer  equation  must  be  treated  as 
an  equivalent  result,  but  it  can  never  be  related  directly  back  to  the  original 
problem,  because  we  have  removed  any  distinctions  between  forward  scattered 
and  unscattered  energy.  Beyond  this  distinction,  the  new  problem  conserves 
energy  in  the  same  degree  as  the  original  problem. 


A. 3  Least  Squares  Derivation 


In  this  section,  we  consider  alternate  methods  for  evaluating  /  and  ci  based  on 
a  least  squares  approach.  We  can  then  compare  results  obtained  using  this  new 
approach  with  the  /  =  Xi_|_i  approach  proposed  by  Wiscombe  (1977).  Let  us 
begin  by  defining  a  least  squares  metric. 


1 

=2tt  J  {P{ii)  -  Ps(/i)}"  diJ.. 

-1 


(.4.13) 
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Using  this  form  of  a  least  squares  equation  requires  that  all  arguments  of  P£;(/i) 
be  square  integrable.  Obviously,  the  S  component  is  not,  so  it  will  need  to  be 


replaced  by  a  function  which  operates  nearly  like  a  delta  function.  That  is,  with 
respect  to  the  component  Legendre  polynomials,  it  will  have  to  act  with  sifting 
properties  almost  like  a  S  function.  As  a  simple  mathematical  device,  in  this 
analysis,  the  S  is  approximated  by  a  Gaussian  form 

(5(/i  -  l)/27r  sa  exp(-6'V2^»^),  (^-14) 

where 

TT 

=  Q~‘^  J  exp{—0^/2g^)sm{0)d0.  (A. 15) 

0 

N(g)  is  a  normalization  term  accounting  for  the  fact  that,  while  we  expect  the 
angular  width  of  the  distribution  about  the  zero  peak  (g)  to  be  narrow,  we  are 
nonetheless  integrating  over  a  unit  sphere  (not  over  a  planar  surface)  such  that 
N(g)  must  be  slightly  greater  than  unity  to  compensate.  By  this  approximation 
of  a  (5  phase  function,  we  are  saying  that  its  width  is  small  with  respect  to 
the  width  of  the  Legendre  components  in  the  expansion.  In  order  to  ensure 
this  feature,  we  restrict  g  to  values  below  some  maximum,  i?,  which  we  have 
arbitrarily  set  at  0.15  radians. 

Using  this  approximation  for  the  S  peak,  can  be  written  as  a  series  of 
six  terms:  C'p  which  represent  the  squaring  of  the  original  phase 

function,  the  S  term,  and  the  Legendre  terms  with  themselves,  and  cross-terms 
between  these  3  basic  types  of  terms.  Of  course,  the  Legendre  terms  of  different 
indices  are  orthogonal  and  integrate  to  zero.  After  some  manipulation,  we  obtain 
the  following  results: 


-1 


where 


£=0 

L 

u< 

i=0 

7r/e 


■2fN{g)Qp{g), 

(.4.16) 

fN\g)Qs{g), 

(.4.171 

L 

=  (1  -  ff 

t=0 

(.4.181 

#  9/\T-.r  /  \n  SlTii  QU^ 

Qi{q)  =  /  exp(— m^/2)  P^[cos(^)m)] - -du^ 


(A.19) 


7r/e 

^  /  X  f  /  9  ,  .  M  sinfpu)  , 

Qp{q)=  /  exp(— M  /2)P[cos(^)m)]  - -du^ 


(A.20) 
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1 


(A21) 


Qsig) 


27rg‘^ 


/exp(-.,^) 

0 


sin(^)M) 

- du. 

Q 


Of  these,  terms  (7o,  Os,  and  represent  the  squaring  of  the  original 
phase  function,  the  Gaussian  peak,  and  the  Legendre  terms  with  themselves, 
respectively.  Ci  represents  the  correlation  between  the  Gaussian  peak  and  the 
original  phase  function,  C2  represents  the  correlation  between  the  Legendre 
terms  and  the  original  phase  function,  and  G4  represents  the  correlation  between 
the  Gaussian  peak  and  the  Legendre  terms.  Because  q  is  small  compared  to  vr, 
it  will  rarely  be  necessary  to  carry  the  Q  integrations  out  over  their  complete 
range. 

is  minimized  by  taking  its  derivative  with  respect  to  the  individual  variables. 
Setting  dZ^  jdci  =  0,  one  obtains  results  from  the  derivatives  of  the  sum  of 
components  G2,  G4  and  G5.  Dividing  by  the  common  factor  2(1  —  f)ai  and 
solving  for  c^,  one  obtains. 


ct 


Xe-fN{g)Qe{g) 

1-/ 


(A.22) 


Notice  that  Qo(g)  =  N  ^{g),  so  as  expected  the  cq  term,  which  carries  all  the 
energy  continues  to  evaluate  to  unity. 

In  the  limit  as  0,  iV  and  the  Qi  both  approach  unity  and  the  same  equation 
is  obtained  for  ci  as  was  obtained  by  the  (5-M  method  [Eq.  (A. 5)].  This  confirms 
that  a  basic  least  squares  approach  is  similar  to  the  (5-M  method.  But,  this 
method  also  provides  information  about  /  and  g^  while  (5-M  only  provided  an 
empirical  formula  for  /. 

We  thus  continue  by  substituting  the  expression  for  Ci  in  the  equation  for  Z^  to 
yield. 


=  Go  -  2/  N{g)  Qp{g)  +  fN\g)  Qs{g)  -  a,  {X,  -  /  N{g)  Qi{g)}\ 

1=0 

(A.23) 

Taking  dZ^/df  =  0  yields. 


fN{g)  = 


r,  /  X  Qp(q)  -  S  aiXiQi(g) 

^n[Q)  _  _ ^=0 _ 

Qsig)  -  aiQ^iig) 

i=0 


(A.24) 


where  Sn  and  Sd  are  numerator  and  denominator  quantities,  respectively. 
Introducing  this  new  form  once  again  to  the  Z^  equation,  one  obtains. 


Z^ 


Co  -  Ko 


sljg) 

Sdig)  ’ 


L 

Ko  =  oiiXf  , 
1=0 


(A.25) 
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where  we  note  that  both  Cq  and  Kq  are  constants  with  respect  to  the  phase 
function  parameter  g.  can  thus  be  minimized  by  varying  g  alone. 

Note  that  we  require  S'^(O)  =  P(l)  —  >  0  for  the  method  to  predict 

any  forward  peak  at  all.  Otherwise,  /  is  set  to  zero  and  ci  =  Xi.  For  phase 
functions  which  pass  the  first  criterion,  g  can  be  determined,  followed  by  /,  and 
finally  the  c^’s. 

To  see  the  properties  of  a  typical  phase  function  when  characterized  by  the 
parameters  /  and  g^  we  evaluated  the  test  particle  size  distribution  at  different 
wavelengths  and  have  plotted  various  functions  of  the  results  combined  with  the 
mode  radius  property  of  the  size  distribution,  the  single  scattering  albedo  ru, 
and  the  wavelength  A,  in  figure  A.l. 

As  expected,  we  found  that  the  ratio  of  g  to  \/ro  is  roughly  constant,  indicating 
that  the  forward  peak  represents  a  diffraction  effect.  And  second,  the  fraction 
/  of  energy  in  the  forward  peak  is  related  to  ru,  because  ru  is  only  reduced  due 
to  energy  absorbed  within  the  particles.  As  more  energy  is  absorbed  inside,  ru 
decreases.  But,  at  the  same  time,  /  increases  as  the  amount  of  diffracted  energy 
remains  the  same,  but  increases  as  a  proportion  of  the  total  scattering. 
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Figure  A.l.  Comparison  of  forward  scattering  parameters  g  and  /  and 
propagation  parameters  ru,  ro  and  A. 


This  analysis  gave  us  considerable  confidence  in  the  least  squares  approach  in 
that  (1)  it  was  able  to  characterize  properties  of  the  forward  peak  region  that 
are  consistent  with  theory  of  scattering,  and  (2)  that  it  produced  calculation 
equations  for  the  ci  that  correspond  to  previous  theory.  However,  there  is  a 


113 


deficiency  in  this  approach  due  to  its  failure  to  produce  positive  definite  phase 
functions  for  low  order  expansions.  A  sample  case  is  shown  in  figure  A. 2,  in 
which  the  predicted  least  squares  curve  for  L  =  2  resulting  from  an  optimal 
/  =  0.54,  produced  a  wide  negative  phase  function  estimate  region  from  7r/2  to 
Svr/d.  The  need  for  a  nonnegative  requirement  was  shown  with  respect  to  the 
scattering  results  presented  in  the  Monte  Carlo  chapter  of  the  main  text.  We 
must  still  look  for  some  representative  value  of  /  which  avoids  negative  values. 


Figure  A. 2.  Results  of  preliminary  least  squares  approach  produce  negative 
phase  function  estimates  at  scattering  angles  between  7r/2  and  Svr/d. 


We  chose  to  work  with  the  L  =  2  (8-stream)  case,  since  it  can  be  solved  in  a 
closed  form.  The  Legendre  expansion  for  this  case  is 

=  (47r)-i [1  +  3ci/i  +  5c2(3/i2  -  l)/2],  (A.26) 


which  is  a  quadratic  equation  in  the  scattering  cosine  /i.  We  can  thus  solve  for 
values  of  /i  which  produce  P'{p)  =  0,  as 


-  _-£L  ^  _£L 

5c2  ^  V  3  ^  25c2 


(A.27) 
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The  requirement  to  avoid  phase  function  zeros  is: 

3c2  +  (5c2  -  1)"  <  1. 


(A.28) 


If  we  then  introduce  the  (5-M  form  for  the  c^,  we  find  that  valid  /  values  will  he 
in  a  range  A/  about  a  mean  value  /,  defined  as, 

-  3Ai  +  2OA2  -  5 
^  ~  18 

^  [25  +  15X1(8X2  -  3X1  -  2)  -  10X2(5X2  + 

^  ~  18 
The  sample  fog  phase  function  data  were  analyzed  on  the  basis  of  Eqs.  (A. 29) 
and  (A. 30)  and  plotted  in  figure  A. 3.  This  figure  reveals  a  section  around  10  /xm 
where  only  the  mean  value  is  plotted  because  A/  has  become  imaginary  over 
this  range.  This  indicates  that  there  is  no  valid  range  of  /  about  the  mean  for 
which  the  phase  function  remains  nonnegative. 


(A.30) 


(A.29) 


Figure  A. 3.  Evaluated  range  of  /  values  (maximum,  minimum,  and  mean)  for 
which  phase  function  remains  nonnegative. 


We  are  thus  led  to  the  conclusion  that  the  traditional  equation  for  evaluating  the 
c^’s  was  insufficient  to  maintain  a  positive  definite  status.  Since  this  equation  is 
based  on  a  standard  least  squares  analysis,  doubt  is  thus  placed  on  this  method. 
An  analysis  of  this  method  and  the  nature  of  standard  phase  functions  reveals 
where  the  problem  arises:  the  least  squares  is  attempting  to  match  the  value  of 
a  function  which  exhibits  a  very  high  peak  at  /i  =  1  and  very  small  values  over 
the  rest  of  its  range.  But,  since  the  least  squares  approach  merely  attempts  to 
minimize  the  distance  of  the  estimate  to  the  curve,  the  large  values  in  the  peak 
region  will  receive  maximum  attention,  and  there  is  no  penalty  for  ‘crossing  the 
line’  into  negative  phase  function  values  at  large  scattering  angles. 
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To  solve  this  situation  we  can  use  a  LLS  optimization: 


1 

=  27r  J  dfj,  {in  [P(jj,)]  —  In  [Pe(jj)]}^  ^  (A. 31) 

-1 


where 


lA(a;) 


ln(a;), 

ln(P^*„)  +  10®(x  -  Pmin), 


®  ^  Pmin 
X  "P:  Pmin 


(A.32) 


where  Pmin  is  chosen  at  some  level  below  the  minimum  value  of  the  tabulated 
phase  function  data.  For  example,  Pmin  =  min(P(/i))/10.  Since  most  phase 
functions  have  minima  in  the  range  a  Pmin  around  10”"^  should  be 

sufficient  to  result  in  only  the  natural  log  function  being  accessed  near  the  end 
of  the  optimization  process,  but  avoid  evaluating  negative  arguments  during  the 
intermediate  stages. 

To  be  compatible  with  the  DOM  we  must  still  describe  our  phase  function 
estimate  using  Eq.  A. 3  since  using  a  LLS  approach  leads  to  the  complications. 
But  then,  it  becomes  impossible  to  analyze  the  least  squares  equation  to 
produce  analytical  solutions  for  the  different  coefficients.  Without  these 
analytical  equations  we  end  up  with  a  multidimensional  search  over  the  ci 
which  could  produce  local  minima  making  the  search  process  very  difficult  and 
time  consuming.  Considering  this  complication,  the  Henyey-Greenstein  phase 
function  (HGPF)  model  was  considered  as  a  model  in  guiding  the  development  of 
an  empirical  method  of  determining  the  ci  values  while  simultaneously  imposing 
a  positive  definite  conditioning  through  the  use  of  the  LLS  equation. 


A. 4  The  Henyey-Greenstein  Phase  Function 


The  Henyey-Greenstein  model  for  the  phase  function  is  well  known  and  has  been 
adopted  by  several  authors  (LeNoble  (1985),  van  de  Hulst  (1980),  Liou  (1980)) 
who  view  it  as  a  model  of  the  typical  behavior  of  aerosols,  in  spite  of  the  fact  that 
the  HGPF  has  no  backscatter  peak  when  the  asymmetry  parameter  is  positive. 

The  use  of  HGPF  is  widespread  because  it  has  a  relatively  simple  analytical 
form  and  because  its  expansion  in  terms  of  Legendre  polynomials  is  particularly 
simple.  The  HGPF  is  written  under  our  normalization  scheme  as. 


PhcAp) 


1  (1-g^) 

dvr  {1  -\-  g'^  —  2gfiyP  ’ 


(A.33) 


where  g  is  the  asymmetry  parameter,  often  set  to  the  value  of  Xi,  and  varies 
between  -1  and  1.  Legendre  coefficients  of  this  distribution  are 

Xhg,£  =  g^-  (A. 34) 


Thus,  ni  g  =  1  the  HGPF  is  a  (5  phase  function. 
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Figure  A. 4.  Estimates  of  the  g  parameter  based  on  the  Henyey-Greenstein 
model.  Xi  for  t'  =  1,  3,  5,  7  were  used  to  estimate  g  for  the  sample  fog  aerosol 
type. 


Most  aerosols  exhibit  Xi  values  that  can  be  approximated  by  some  g  value 
between  0.85  and  1.00,  as  illustrated  in  figure  A. 4.  There,  we  used  the  Xi  to 
compute  g  estimates  using  the  equation  g  sa  X^^ .  The  resulting  estimates  seen 
in  this  figure  show  some  dispersion  yet  reveal  that  the  HGPF  is  a  realistic  model. 

This  resemblance  between  the  HGPF  and  aerosol  phase  function  coefficients 
is  useful  in  providing  some  insights  into  problems  encountered  when  finite 
expansions  of  this  form  are  needed.  We  identify  this  truncated  form  using  the 
function  Phg{jJ^)^ 


L 

Phg{p)  ~  Phg{p)  =  ^  o:£  g^  Pt{}i).  (A. 35) 

The  family  of  truncated  HGPF  functions  behaves  very  differently  from  the  full 
HGPF.  For  example,  figure  A. 5  shows  the  behaviors  of  several  approximations 
to  a  delta  function  for  different  even  orders  of  expansion  {L  =  2, 4, 6, 8). 
From  the  figure  we  see  a  ringing  behavior  similar  to  the  Gibbs  phenomenon 
observed  in  Fourier  series  expansions  around  functional  discontinuities.  These 
are  particularly  severe  when  modeling  a  8  function,  and  since  most  aerosols 
exhibit  highly  forward  scattering  behaviors  similar  to  (5’s,  we  should  expect 
ringing  there  as  well. 

To  remove  these  artifacts  and  yet  maintain  the  HGPF  form,  we  find  that  we 
must  reduce  the  g  value  down  to  some  threshold  maximum  value  that  cannot  be 
exceeded  without  producing  a  negative  phase  function.  Figure  A. 6  shows  plots 
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Figure  A. 5.  Legendre  expansion  approximations  to  the  S  phase  function  at 
different  orders  of  expansion,  showing  zero  crossings  for  orders  L  =  2,  4,  6,  and 

8. 

of  a  few  of  these  truncated  HGPF  forms  at  maximum  g  values.  We  have  plotted 
these  maximum  g  values  out  in  figure  A. 7.  In  this  same  figure,  we  have  also 
plotted  out  the  results  of  a  second  investigation  where  we  asked  the  question: 
What  if  we  perform  a  (5-M  correction  on  a  truncated  HGPF  representation?  This 
correction  involves  setting  /  =  leading  to, 

_  / 

1  —  g^+^ 

This  second  curve  indicates  that  when  we  impose  a  (5-M  forward  peak  truncation, 
we  increase  the  effective  maximum  g  that  can  be  treated  within  the  model. 
But,  comparing  these  new  maxima  with  characteristic  g  values  of  a  real  aerosol 
(figure  A. 4)  reveals  that  we  cannot  adequately  characterize  these  aerosols  at  L 
values  less  than  about  16  or  18  using  the  (5-M. 

Since  using  an  T  of  16  or  18  would  impose  heavy  memory  burdens  on  the  RT 
computations,  some  form  for  the  ci  computation  was  needed  that  provided  a 
system  similar  to  that  originally  developed,  yet  avoided  the  negative  phase 
function  problem.  As  an  initial  solution,  it  was  observed  that  if  one  started 
with  an  HGPF  characterized  by  a  value  g,  then  this  value  must  be  corrected 
downward  to  a  value  gi^  corresponding  to  the  value  of  g  for  order  L  which  is 
the  maximum  g  parameter  that  produces  a  positive  definite  expansion  at  that 
order.  To  obtain  this  reduction,  one  would  need  to  multiply  expansion  coefficient 
7  by  a  factor  e^,  where  e  <  1.  But  if  we  started  with  the  expansion  coefficients 
Cl  =  {Xi  —  /)/(!  —  /),  these  coefficients  could  be  viewed  as  a  series  of  values 
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Figure  A. 6.  Legendre  expansion  approximations  to  the  maximum  positive- 
definite  HGPF’s  at  different  orders  of  expansion.  Amplified  region  shows  curves 
approach  zero  due  to  local  minima  closest  to  /i  =  —1. 


1.0 


OJ 


0.9 


03 

PLh 

>5 

S-l 


0.8 


0.7 


>5 

CO 


0.6 


X! 

03 


0.5 


0.4 


o 

For  6-M  Mod.  to  HG  Di.st. 

O 

For  Truncated  HG  Distrib. 

. o 

..-■O . ° 

o . ®  ■■ 

..0 

.....  O 

..O  ' 

...O'' 

o  ■ " 

0  2  4  6  8  10  12 

Legendre  Summation  Limit 


Figure  A. 7.  Limiting  maximum  value  allowable  for  g  when  imposing  a  positive 
definite  condition  on  an  order  L  Legendre  approximation  of  the  Henyey- 
Greenstein  phase  function,  and  when  using  the  (5-M  correction  method. 

which  decrease  like  a  Henyey-Greenstein  set,  which  leads  to  the  proposition  of 
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a  new  empirical  relationship  for  the  ci  of 


[X,-fN{g)Q,{g)]  , 

[1  -  /] 


(A.36) 


Here  we  have  used  the  standard  least  squares  results  and  augmented  them  by 
the  appended  factor.  Since  Xq  =  1  and  N { g)  Q g)  =  1,  we  still  have  cq  =  1, 
which  conserves  scattered  energy. 

To  understand  this  solution  somewhat  better,  consider  the  approximation  that 
the  (5-M  is  making:  the  goal  under  the  (5-M  is  to  exactly  match  the  Legendre 
expansion  coefficients  from  the  original  series  expansion  in  the  new  expansion, 


X,  =  +  f)c,. 


(A.37) 


Then,  by  using  values  for  (5^  =  1,  the  standard  equation  is  obtained  for  c^.  In  our 
method  we  have  Eq.  (36)  for  the  c^’s  and  wish  to  work  backwards  to  compute 
the  coefficients  8f. 


X,  =  f  8, +  e^  [X,  -fN{g)  Qi{g)i  (A.38) 

which  leads  to, 

8i  =  y  (1  -  +  N{g)Qi{g)  e^.  (A.39) 

Obviously,  Eq.  (39)  produces  81  =  1  whenever  e  =  1  and  =  0.  And  these 
two  parameters  approach  these  limiting  values  as  the  order  L  of  the  expansion 
increases.  Thus,  the  (5-M  technique  can  be  seen  to  model  the  8  function  using 
a  truncated  expansion  of  the  exact  8  function  coefficients.  Our  method  models 
the  8  using  optimized  coefficients  based  on  the  order  of  the  expansion.  The 
difference  is  that  the  current  solution  uses  coefficients  that  minimize  the  ringing 
behavior  in  the  truncated  series  results. 

We  show  the  results  of  the  LLS  analysis  in  figures  A. 8  and  A. 9.  In  figure  A. 8, 
the  optimal  value  of  e  and  the  resulting  are  plotted  for  an  T  =  2  case.  Note 
that  e  is  a  maximum  when  the  sum  of  squares  is  a  minimum,  indicating  that 
there  appears  to  be  some  sense  in  which  the  parameterization  of  e  produces 
maximum  values  near  the  solution  point,  leading  us  to  conclude  that  (1)  the 
method  is  stable,  and  (2)  the  method  attempts  only  to  remove  the  negative 
effects  rather  than  abnormally  reducing  the  influence  of  higher  Legendre  terms. 

In  figure  A. 9,  the  expansion  resulting  from  this  choice  of  e-f  pair  is  plotted  for 
L  =  2  and  for  a  set  of  results  for  L  =  4.  In  the  L  =  4  case,  /  =  0.58  and  e  =  0.93. 
As  seen,  the  L  =  4  expansion  captures  the  forward  hemisphere  behavior  better 
than  the  L  =  2  expansion. 
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Forward  Scattering  Fraction 

Figure  A. 8.  The  minimum  error  value  of  e  is  plotted  along  with  the  scaled 
squared  error  associated  with  that  value  as  a  function  of  the  forward  scatter 
fraction  /. 


Figure  A. 9.  Comparison  of  Legendre  expansions  for  L  =  2  and  L  =  4  using 
positive  definite  choices  for  /  and  e  against  calculated  phase  function  data. 


121 


A. 5  Conclusions 

From  figure  A. 3  it  was  shown  that  a  simple  rule  of  thumb  for  modifying  the  ci 
coefficients,  based  on  the  forward  scattering  fraction  /  alone,  can  lead  to  negative 
phase  function  evaluations  at  some  angles.  We  thus  introduced  an  empirical 
equation  for  C£  that  employed  the  modification  by  the  factor  e^.  This  new  form 
attenuates  higher  order  Legendre  coefficients  to  ensure  a  positive-definite  phase 
function  and  remove  the  forward  peak. 

The  method  developed  here  results  in  a  much  more  robust  means  of  obtaining  an 
optimal  approximate  phase  function  including  the  forward  scattering  correction, 
while  simultaneously  avoiding  the  problems  inherent  with  negative  phase 
function  predictions.  Chapter  3  illustrated  that  these  new  techniques  produce 
better  flux  predictions  than  the  (5-M  technique  of  Wiscombe. 

In  treating  scattering  problems,  once  an  order  of  expansion  has  been  chosen 
for  a  simulation,  the  phase  function  properties  can  be  preprocessed  to  derive 
appropriate  /,  and  e  values  that  lead  to  modifications  to  the  extinction 
coefficient  cr,  the  single  scattering  albedo  ru  and  the  modified  Legendre 
coefficients  ci  that  describe  that  aerosol. 
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STEDP  MT  M 
ATTN  MR  BOWERS 
DUGWAY  UT  84022-5000 

PL  WE  1 

KIRTLAND  AFB  NM  87118-6008 

USAF  ROME  LAB  TECH  1 

CORRIDOR  W  STE  262  RL  SUL 
26  ELECTR  PKWY  BLD  106 
GRIFFISS  AFB  NY  13441-4514 

AFMC  DOW  1 


WRIGHT  PATTERSON  AFB  OH  45433-5000 


126 


ARMY  FIELD  ARTILLERY  SCHOOL  1 

ATSF  TSM  TA 
FT  SILL  OK  73503-5600 

NAVAL  AIR  DEV  CTR  1 

CODE  5012 
ATTN  AL  SALIK 
WARMINISTER  PA  18974 

ARMY  FOREIGN  SCI  TECH  CTR  1 

CM 

220  7TH  STREET  NE 
CHARLOTTESVILLE  VA  22448-5000 

NAVAL  SURFACE  WEAPONS  CTR  1 

CODE  G63 

DAHLGREN  VA  22448-5000 

ARMY  OEC  1 


CSTE  EFS 
PARK  CENTER  IV 
4501  FORD  AVE 
ALEXANDRIA  VA  22302-1458 


ARMY  CORPS  OF  ENGRS  1 

ENGR  TOPOGRAPHICS  LAB 
ETL  GS  LB 

FT  BELVOIR  VA  22060 

ARMY  TOPO  ENGR  CTR  1 

CETEC  ZC  1 

FT  BELVOIR  VA  22060-5546 

SCIENCE  AND  TECHNOLOGY  CORP  1 

101  RESEARCH  DRIVE 
HAMPTON  VA  23666-1340 

ARMY  NUCLEAR  CML  AGCY  1 

MONA  ZB  BLDG  2073 
SPRINGFIELD  VA  22150-3198 

USATRADOC  1 

ATCD  FA 

FT  MONROE  VA  23651-5170 

ARMY  TRADOC  ANALYSIS  CTR  1 

ATRC  WSS  R 
WSMR  NM  88002-5502 

ARMY  RESEARCH  LABORATORY  1 

AMSRL  IS  EW 
BATTLEFIELD  ENVIR  DIV 
WSMR  NM  88002-5501 

ARMY  RESEARCH  LABORATORY  1 

AMSRL  IS  ES 

BATTLEFIELD  ENVIR  DIV 
ADELPHI  MD  20783-1145 
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DTIC 

8725  JOHN  J  KINGMAN  RD 
STE  0944 

FT  BELVOIR  VA  22060-6218 

ARMY  MISSILE  CMND 
AMSMI 

REDSTONE  ARSENAL  AL  35898-5243 

ARMY  DUGWAY  PROVING  GRD 
STEDP3 

DUGWAY  UT  84022-5000 

USTRADOC 
ATCD  FA 

FT  MONROE  VA  23651-5170 

WSMR  TECH  LIBRARY  BR 
STEWS  IM  IT 
WSMR  NM  88002 

US  MILITARY  ACADEMY 
MATHEMATICAL  SCI  CTR  EXCELLENCE 
DEPT  OF  MATHEMATICAL  SCIENCES 
ATTN  MDN  A  (MAJ  DON  ENGEN) 
THAYER  HALL 
WEST  POINT  NY  10996-1786 

ARMY  RESEARCH  LABORATORY 
ATTN  D  TOFSTED 
AMSRL  IS  EW 
BATTLEFIELD  ENVIR  DIV 
WSMR  NM  88002-5501 


Record  Copy 


TOTAL 


